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Возможности оценки динамических состояний железнодорожных [n] [n] 
I 
транспортных средств: структурное математическое моделирование : 
P.C. Большаков ‚В.Е. Гозбенко@©, К.Ч. Выонг© [=] 
Иркутский государственный университет путей сообщения, г. Иркутск, Российская Федерация 
bolshakov_rs@mail.ru EDN: PEUPEC 
Аннотация 


Введение. Увеличение скоростей движения железнодорожного транспорта и повышение нагрузок на оси колес- 
ных пар обуславливают необходимость модернизации существующего парка. Научные исследования в области 
динамики подвижного состава направлены на учёт колебательных процессов, возникающих при движении же- 
лезнодорожных транспортных средств в традиционном конструктивном исполнении. Присоединение дополни- 
тельных элементов рассматривалось на уровне сцепки двух вагонов и присоединении третьей тележки в центре 
тяжести железнодорожной платформы. Построению математических моделей, позволяющих оценить динамиче- 
ские состояния таких конструктивных решений, в научной литературе не уделено достаточно внимания. Цель 
данного исследования — создать метод оценки динамических состояний вагона. Рассматривается ситуация, ко- 
гда в его структуру вводится дополнительная совокупность масс-инерционных и упругих элементов, причем от 
корректировки их параметров зависит общее динамическое состояние транспортного средства. 

Материалы и методы. Базовым инструментом проведения исследований является структурное математическое 
моделирование, в основе которого лежит подход, когда исходная расчетная схема представляет собой механиче- 
скую колебательную систему в виде твердого тела на упругих опорах с дополнительной введёнными в её струк- 
туру типовыми элементами. Динамическим аналогом используемой расчетной схемы является структурная схема 
системы автоматического управления, применение которой позволяет детализировать связи между типовыми 
упругими и масс-инерционными элементами. 

Результаты исследования. Предложен метод оценки динамических состояний железнодорожных транспортных 
средств, основанный на построении математических моделей, с учетом введения дополнительной структуры 
масс-инерционных и упругих элементов. Исследовано влияние дополнительных параметров на динамическое 
состояние транспортного средства. Получены аналитические соотношения, позволяющие при изменении соот- 
ветствующих параметров технического объекта снизить динамические нагрузки на основные конструктивные 
упругие элементы. Приведена передаточная функция межпарциальных связей, позволяющая контролировать вза- 
имодействие между координатами движения транспортного средства при действии двух кинематических возму- 
щений синфазного типа. 

Обсуждение и заключение. Сформированная математическая модель позволяет оценить динамическое состоя- 
ние железнодорожного транспортного средства в условиях действия кинематических возмущений. Результаты 
исследований могут быть применены при модернизации существующих и создании новых транспортных средств 
с улучшенной динамикой. 


Ключевые слова: динамика транспортных средств, динамическое состояние, структурное математическое 
моделирование, структурная схема 
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Assessment of Dynamic States of Railway Vehicles: Structural Mathematical Modeling 
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Abstract 

Introduction. The speed rise of railway transport and an increase in the loads on the axles of wheelsets necessitate the 
modernization of the existing fleet. Scientific studies in the field of rolling stock dynamics are aimed at taking into account 
the oscillatory processes that occur during the movement of railway vehicles in a traditional design. The attachment of 
supplementary elements was considered at the coupling level of two cars and the attachment of a third trolley in the center 
of gravity of the railway platform. The scientific literature has not paid enough attention to the construction of 
mathematical models that make it possible to assess the dynamic states of such constructive solutions. The objective of 
this study is to create a method for evaluating the dynamic conditions of a car. The situation is considered when an 
additional set of mass-inertial and elastic elements is introduced into its structure, and the general dynamic condition of 
the vehicle depends on the adjustment of their parameters. 

Materials and Methods. The basic research tool is the structural mathematical modeling, which is based on an approach 
where the source design scheme is a mechanical oscillatory system in the form of a solid body on elastic supports with 
supplementary typical elements introduced into its structure. The dynamic analogue of the calculation scheme used is the 
block diagram ofthe automatic control system, the use of which provides detailing the connections between typical elastic 
and mass-inertia elements. 

Results. ^ method for estimating the dynamic states of railway vehicles is proposed. It is based on the construction of 
mathematical models, taking into account the introduction of an additional structure of mass-inertia and elastic elements. 
The impact of additional parameters on the dynamic condition of the vehicle is investigated. Analytical relations have 
been obtained that provide reducing the dynamic loads on the major structural elastic elements when changing the 
corresponding parameters of a technical object. The transfer function of interpartial relations 1s given, which provides 
controlling the interaction between the coordinates of the vehicle movement under the action of two kinematic 
disturbances of the in-phase type. 

Discussion and Conclusion. The generated mathematical model provides for assessment, monitoring and control of the 
dynamic state of the vehicle under conditions of kinematic disturbances. The research results can be used to modernize 
existing vehicles and create new ones with improved dynamics. 


Keywords: vehicle dynamics, dynamic condition, structural mathematical modeling, block diagram 
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Введение. Расширение технико-экономических связей межрегионального уровня, обеспечение роста произ- 
водственно-технического потенциала страны и поддержание разрастающейся системы международных торгово- 
коммерческих отношений во многом зависит и опирается на железнодорожный транспорт [1]. Объёмы перевозок, 
выполняемых железнодорожным транспортом, постоянно возрастают, что обусловливает необходимость учета 
нежелательных эффектов от повышения динамических нагрузок [2], напрямую влияющих на надежность эксплу- 
атации как подвижного состава, так и верхнего строения пути. Несмотря на наличие таких негативных факторов, 
необходимо выполнять плановые показатели, к которым можно отнести повышение участковой скорости, со- 
блюдение норм веса и длины составов, увеличение осевой нагрузки до 30 тонн и более. Это отражает реальные 
запросы развития российской экономики и стимулирует создание новых более мощных локомотивов, обновление 
парка подвижных средств, модернизацию путевого хозяйства [3]. Вместе с тем нельзя не рассматривать и BO3- 
можности проявления негативных последствий интенсификации перевозочных процессов. Одним из важных во- 
просов является возрастание или увеличение темпов износа верхнего строения пути с соответствующими выте- 
кающими сложностями [4]. 
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В текущий момент времени повышенное внимание уделяется развитию методологии оценки динамического со- 
стояния подвижного состава, взаимодействию технических средств с рельсовыми путями, экономии электроэнер- 
гии, повышению надёжности и безопасности перевозочных процессов. Методология математического моделирова- 
ния описана, например, в [5]. Вместе с тем, существуют и другие возможности для поиска рациональных реше- 
ний [6]. Важно уделить внимание вопросам модернизации существующего парка грузовых вагонов, эксплуатация 
которых при повышенных нагрузках уже не является эффективной [7]. Одним из подходов, который мог бы быть 
принят к разработке, является концепция установки для грузовых 4-осных вагонов дополнительной двухосной те- 
лежки [8]. В этом случае можно ожидать более равномерного распределения нагрузки на верхнее строение пути, а 
также возможность увеличения веса перевозимых грузов при сохранении регламента на значения осевой нагрузки 
в пределах 22 т [9]. Неравномерность параметров контакта колесо-рельс инициирует колебательные движения Ba- 
гона, что формирует в свою очередь колебательные движения корпуса вагона. Процесс колебаний формируется 
также условиями взаимодействия вагонов внутри состава [10]. В этом случае возникают динамические реакции 
связей, накладывающиеся на статические составляющие полной реакции, что может существенно увеличить уро- 
вень динамических взаимодействий в контакте колесо-рельс [11, 12]. Однако возможностям структурного матема- 
тического моделирования при оценке динамических состояний железнодорожных транспортных средств при вве- 
дении дополнительных связей пока не уделялось достаточно внимания. Поэтому целью данного исследования яви- 
лось формирование метода оценки динамических состояний железнодорожного транспортного средства при введе- 
нии в его структуру дополнительной совокупности масс-инерционных и упругих элементов, корректировка пара- 
метров которых будет влиять на общее динамическое состояние транспортного средства. 

Материалы и методы. Методологической основой исследований выбрана структурная теория виброзащит- 
ных систем, позволяющая в линейной постановке с учетом сосредоточенных параметров и малых колебаний от- 
носительно положения статического равновесия или установившегося процесса, достаточно точно оценить ди- 
намические свойства железнодорожного транспортного средства. Расчётной схемой является механическая ко- 
лебательная система с динамическим эквивалентом в виде структурной схемы системы автоматического управ- 
ления. Это позволяет детализировать связи между элементами транспортного средства, а также использовать 
методы, характерные для теории автоматического управления (передаточные функции, преобразования струк- 
турных схем, свертки и упрощения, частотные характеристики) [13]. 

Стандартные конструктивные схемы грузовых железнодорожных транспортных средств описываются извест- 
ными расчётными схемами, а их динамические особенности можно оценить с использованием линейных расчёт- 
ных схем. Рассматривается железнодорожное транспортное средство в виде грузового четырехосного вагона, 
предназначенного для перевозки тяжёлых грузов, таких как каменный уголь, руда, песок, металлопрокат, мало- 
габаритные металлические конструкции и др. (рис. 1). 
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Рис. 1. Принципиальная схема грузового четырехосного вагона 


Структура представленного железнодорожного транспортного средства содержит кузов, обладающий мас- 
сой М, моментом инерции Ји опирающийся на две двухосные тележки, условно представленные как совокупно- 
сти масс-инерционных и упругих элементов. Динамические особенности рассматриваемой системы показывают 
наличие избыточного влияния на конструктивные элементы тележек железнодорожного транспортного средства. 
В этой связи для улучшения его динамических СВОЙСТВ В структуру подвески дополнительно вводится третья 
двухосная тележка, которая размещается в центре вагона (рис. 2). 
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Рис. 2. Принципиальная схема грузового шестиосного вагона 


При собственном весе вагона около 23 тонн и осевой нагрузке в 22 тонны, это может обеспечить увеличение 
веса перевозимых грузов ориентировочно на 20 тонн, то есть существенно повысить эффективность перевозоч- 
ных процессов, не создавая излишних динамических нагрузок на верхнее строение пути (ВСП). Практическая 
реализация предлагаемого подхода требует модернизации конструкции типового вагона, которая заключается в 
создании узла крепления дополнительной тележки для уменьшения нагрузки на ось грузового вагона за счёт 
введения дополнительной двухосной тележки. Модернизация узла крепления осуществляется в тех же конструк- 
тивно-технических формах, что и крепления с использованием шкворней в двух «штатных» двухосных тележках, 
то есть через установку надшкворневой и подшкворневой балок с соответствующим шкворневым узлом стан- 
дартного типа. 

Предлагаемая методика повышения эффективности использования четырехосных грузовых вагонов в усло- 
виях дополнительных нагружений ориентирована на решение проблемы повышения нагрузки на ось колёсной 
пары до 30 тонн и увеличении скорости движения поездов. Это достигается модернизацией типового грузового 
вагона путём установки дополнительной двухосной тележки соответствующим устройством, обеспечивающим 
условия её динамического взаимодействия с конструкцией рамы грузового вагона. 

Установка дополнительной двухосной тележки за счёт перераспределения нагрузки между общим набором 
колёсных пар обеспечивает возможности перевозки грузов больших по весу при одновременном снижении 
нагрузки на ось. Это обеспечивает более рациональные условия эксплуатации рельсового пути и верхнего стро- 
ения железнодорожного полотна при сохранении приемлемой длины состава. 

Предполагается, что шкворневый узел будет иметь упругую резиновую прокладку, обеспечивающую амор- 
тизацию взаимодействия тележки и рамы корпуса вагона. При этом не предполагается модернизация дисков ко- 
лесных пар. Изменяются лишь форма проводки элементов «штатной» пневматической системы торможения со- 
става и конфигурация проводки пневмопроводов. 

Результаты исследования. Расчётная схема рассматриваемого железнодорожного транспортного средства B 
первом приближении может быть представлена в виде механической колебательной системы, состоящей из твёр- 
дого тела, обладающего массой т, и моментом инерции J, опирающемся на три упругих элемента с жесткостями kı, 


ko, k2. Кинематические воздействия представлены синфазными гармоническими функциями одной частоты (рис. 3). 


У M,J 


21 


Рис. 3. Расчётная схема железнодорожного транспортного средства 
при кинематическом возмущении (21, 22, 20) 
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Центр масс системы — точка О находится на расстояниях lı и l2 от концов твёрдого тела (точек 41, В). Эле- 
мент с жесткостью Áo закреплён в точках D и Di; расстояние OD, обозначено как lo. Движение системы рассмат- 
ривается в координатах yi, y» и yo, ф, связанных с неподвижным базисом. В расчётах используются соотношения: 

Yo = ay, +by2, ф = с(у› ул), Yı = yo - hi$; ya = yo +, 
1 l 1 (1) 
2 1 
Ур = У tho, a= b= C= . 
1+1, 1+1, L+ 


Для вывода дифференциальных уравнений движения используется формализм Лагранжа [14], что требует по- 
строения выражений для кинетической и потенциальной энергий. В данном случае имеем: 


2 


T = т(ау кру) tJe? (у-у), (2) 


1 1 1 
П= № (у а)? + 5 ko (Ур E | 5 hb» zy. (3) 


Учтём, что потенциальная энергия (3) может быть записана также в виде: 


1 1 1 
П E (yi -aJ tu (> -2,) +5 [ayı + by; t lhc(y; -a)-a] Е 


1 2 1 32 d] 2 (4) 
Pu (yı -z) rus (уг —22) Ко [в +b yo –2, | > 
где ai =а- loc, bi =b + he. 
Система уравнений движения в координатах у, y» во временной области примет вид: 
я (та? + Је? )+ У (x, +kya?)- y? (Je? -mab) + yokoajb, = kızı + Коа120, (5) 
уз (mb? + Jc”) + ys (ky + kobe )— yi (Jc? —mab) + уа = kzz, + kobizo. (6) 


После применения интегральных преобразований Лапласа при нулевых начальных условиях [15] система 
уравнений (5), (6) может быть представлена в операторной форме: 


У (та? + Jc? )p? +k, + koa |- У» [ (ve? -mab) p? — koayby | = 2 + 00120, (7) 
y; | (mb? +?) p? +k + kyb? | E [ze -mab) р? -koad | = 2, КЬ, (8) 


где р = јо — комплексная переменная ( / = 4/—1 ), значок <-> над переменной означает её изображение по Лапласу [5]. 

На основе (7), (8) построена, в соответствии с положениями метода структурного математического модели- 
рования [5], структурная математическая модель в виде структурной схемы эквивалентной в динамическом от- 
ношении системе автоматического управления (рис. 4). 


(Je – Mab) p? – k,ajb, 


1 1 \ 


2_ 2_ x 
^ (Ма? + Јс?)р? +k, + koa? ао de (Mb? + Jc?) p? + ky + kob? Va 
2 m z 
k «—*! pe — 
k а 4 m L| 20 
061 Кор! = 


Рис. 4. Структурная математическая модель железнодорожного транспортного средства 


Обсуждение и заключение. Особенности системы заключаются в том, что она имеет два парциальных блока, 
каждый из которых определяет соответствующие парциальные частоты: 


2 
2 ki+koai 


ny -—————, 9 

i ma? + Jc? 9) 
ky + kob? 

Э.а. 2 071 

7? > nbi + Je?” (10) 


которые, в свою очередь, определяют границы расположения частот собственных колебаний системы в целом: 
2 2 2,2 2 
Особ < И < 15 < 05, (11) 
ГДЕ W1co6, (02c06 — частоты собственных колебаний системы, при работе на которых возможно возникновение pe- 
зонансных режимов. 
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К числу особенностей системы относится наличие нескольких одновременно действующих внешних возму- 
щений. Полагая в целях упрощения, чтот 2, = 2, = Zo (вполне допустимо на стадиях предварительных динами- 


ческих оценок) примем, что на вход первого и второго парциальных блоков действуют силовые факторы: 
Qi =7(k,+koa,), (12) 
О, =Z(k, +k ob). (13) 
Используя структурную схему Ha рис. 4, запишем выражения для передаточных функций, полагая, что между 
внешними факторами возмущения имеется связь, формируемая соотношением: 


Q = Q, 0, = aQ, (14) 
Риа (kı + koa; )| (mb? Jc?) p? +, + kob? | e o (Ks +в, | (Jc? - mab) р? сһаһ). T 
2 А(р) 
w,(p)- y; _ a(k, + kobi )| (ma? +? р? +k, + koa? |+ (kı + koa, )| (Je? — mab) р? cha | ag 
2 А(р) 
где 
A(p) = [ (ma? +?) p? +k, + kya? || (mb? +?) p? +k, + kyb? |-| (Je? - mab) p? -komb |, (17) 


является частотным характеристическим уравнением системы. 
Числителем передаточных функций в выражениях (15), (16) определяют режимы динамического гашения ко- 
лебаний, что может быть детализировано из уравнений, получаемых при «обнулении» числителей (15), (16): 
2 
р (ki + оа) (К + kob? ) - a (Ks + kobi Кое a 
O |дин = Й 
(ki + а; ) (ть? + Jc? ) +a(k, + kob, )( Jc? - mab) 


"A (kz + kobi ) (в + koa? ) - o (I + koa, )koaib T 
UP a (kika ko (та? + Je?) (ik; + koay (Jc? - mab) 


Из выражений (18), (19) следует, что в системе с двумя степенями свободы при наличии двух связанных 
между собой возмущающих факторов возможно возникновение режимов динамического гашения колебаний на 
двух частотах, параметры которых зависят от значений коэффициента связности а. Этот коэффициент может 
принимать отрицательные, положительные и нулевые значения. 

Из анализа структурной математической модели (рис. 4) следует также и возможность возникновения особого 
динамического режима на частоте: 

Е kab; 

men Je? —mab’ 
когда парциальные блоки получают возможность разъединения. В этом случае система (рис. 3) распадается Ha 
два автономных блока, которые не будут создавать ситуаций взаимодействия парциальных структур. 

Реализация такого режима может привести к существенной разнице в значениях отклонений в точках 41 и В! 
и «разбросу» значений динамических реакций в точках A), Bı и Dj. Определение динамических реакций в точках 
А, B u D может быть реализовано по методике, изложенной в работе [14], в которой динамическая реакция опре- 
деляется как произведение динамического смещения (например, точках 41, Bı и Di) на значение приведённой 


(20) 


динамической жёсткости. 
Для координаты yı, динамическое смещение определяется из выражения для передаточной функции (15), а 


для координаты у,, из выражения (16). Для определения приведённой динамической жёсткости используется 
частотное характеристическое уравнение (17): 


[ (ve? —mab) p? -kyab] 


kaa, 
Kap (p) =k, + — . (21) 
d ( ) i (mb? +?) р? +k, +kob? (mb? +?) р? +k, + kob? 
Аналогично может быть найдено значение приведённой динамической жёсткости по координате у, 
2 
2 2 
Е Б?Ь (Jc -mab) p — kab; 
про (р) = ky F о : | | (22) 


(та? +?) р? +k, + koa? (та? +?) р? +k, + koa? | 
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Для определения конкретных значений k, npl ( p) k ( p) необходимо найти значение частоты à», которая 
определяется из условия, что y; / y, =1. В более общем случае предполагается, что: 
22 1, (23) 
У 
где ү — коэффициент связности движения элементов по координатам у, и yı. Таким образом этот коэффициент 
(для случая y = 1) может быть записан в виде: 
y 9 (k, + kob; )| (ma? + Jc? )р? + + koa? | + (А, + Коа ) (ve? — mab) p- koayby | 
У (k + Коа} )| (me? t Jc? ) 2? +k, + kb? |+ a (Ks + kb, ) (ze? —mab) p? ар | 


Q4) 


Принимая конкретные значения Y, можно найти частоты колебаний для движения рассматриваемого объекта 
в координатах у, и yı. Например, при y = | частота поступательных вертикальных колебаний твёрдого тела опре- 
делится выражением: 
у 95 kb ) (к + koa? ) - (ki + оа) Koaiby + (Ki + koa ) koi — 
n a(k, + kob, (ma? + Jc? )+ (k; + koa; (Jc? - mab) — 
(А + koa; (ks + kob? )- (kz + kobi ) (Ки + koa? ) + а ( + kobi ) а 


-(k + koa, (mb? + Je? ]- a (E; + kobi ) (<? - mab) 


(6) 


(25) 


Частота o синфазных гармонических возмущений 21(0, 22(1) и zo(f) обеспечивает движение твёрдого тела при 


Фф = 0, то есть y, = У, = yp. При других значениях y (у = 1) определяются y, и у,, а Ha их основе из геометриче- 
ских соображений определяется значение yp. Для твёрдого тела при известных y; и у, легко может быть 


найдено положение центра вращения (или колебаний) [12]. 


Динамические реакции связей R4, Кр и Rp, могут быть в первом приближении найдены по формулам: 
Ra = hy; Ки = КУ, Кр = koyp. 
В общем случае динамические реакции определяются с использованием динамических смещений у, У, ур, 


определяемых выражениями (15), (16). Что касается динамического смещения в точке D, то используется выра- 


жение ур =a,y,+b,y,, параметры которого определяются вышеприведёнными значениями. В выраже- 


ния (15), (16) могут быть введены данные о связности силовых факторов воздействия (параметр a). Для получе- 
ния конкретных данных о значениях динамических реакций вводятся параметры частоты, при которой реализу- 


ется соотношение амплитуд колебаний y; и y, через коэффициент связности амплитуд колебаний y. 


Полная реакция связей в точках A, В и D определяется суммой двух компонент: статической и динамической. 
Статическая компонента может быть найдена из выражения для передаточных функций динамических реакций 
при р = 0 и задании параметров статической нагрузки (вес вагона и перевозимого груза). При интенсивном pa3- 


витии колебательных процессов, когда возрастают колебания по координатам yr Уз А Ур, полная реакция может 


изменяться в значительных пределах и отличаться от статической реакции связей. При наличии динамической 
компоненты полная реакция может принимать различные значения, в частности, нулевые или отрицательные. 

Сформированная в рамках предложенного метода математическая модель, обозначенная выражением (25), 
позволяет оценивать динамическое состояние железнодорожных транспортных средств при введении дополни- 
тельных связей в их структуру для формирования комплекса рекомендаций по получению устойчивых режимов 
их эксплуатации. Изучение особенностей системы при помощи подходов, характерных для структурного мате- 
матического моделирования, позволяет детализировано рассмотреть связи между элементами. Применительно к 
рассматриваемому техническому объекту это даёт возможность корректировки динамического состояния техни- 
ческого объекта на основе варьирования параметров совокупности дополнительно введённых элементов и поз- 
воляет уменьшить нагрузку на основные части подвески, а также установить наличие в системе собственных 
частот и частот динамического гашения колебаний. 
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В дальнейшем планируется проведение исследований при введении в структуру транспортного средства 
демпферов и устройств преобразования движения для оценки возможностей структурного математического мо- 
делирования. Также интересным направлением является оценка возможностей изменения динамических реакций 
в зависимости от внешнего воздействия, что позволит оценить усилия, прилагаемые к различным элементами 
подвески транспортного средства. 


Список литературы / References 

1. Коссов B.C., Князев Д.А., Красюков Н.Ф., Махутов H.A., Гаденин М.М. Нормативная база обеспечения 
безопасной эксплуатации железнодорожной техники по ресурсу несущих конструкций. Мир транспорта. 
2023;21(3):106—114. https://doi.org/10.30932/1992-3252-2023-21-3-10 

Kossov VS, Knyazev DA, Krasyukov NF, Makhutov NA, Gadenin MM. Regulatory Framework for Ensuring the 
Safe Operation of Railway Equipment Based on the Service Life of Load-Bearing Structures. World of Transport and 
Transportation. 2023;21(3):254—262. https://doi.org/10.30932/1992-3252-2023-21-3-10 

2. Махутов H.A., Лапидус b.M., Гаденин М.М., Титов Е.Ю. Задачи и перспективы развития научных 
исследований в рамках сотрудничества между ОАО «РЖД» и Российской академией наук. Железнодорожный 


транспорт. 2023;(7):6. 

Makhutov NA, Lapidus ВМ, Gadenin MM, Titov EYu. Tasks and Prospects for the Development of Scientific 
Research within the Framework of Cooperation between JSC “Russian Railways” and the Russian Academy of Sciences. 
Zheleznodorozhnyi transport. 2023;(7):6—11. (In Russ.). 

3. Лапидус Б.М. Задачи опережающего развития российских железных дорог. Железнодорожный транспорт. 
2023;(2):4—14. 

Lapidus BM. Tasks of Advanced Development of Russian Railways. Zheleznodorozhnyi transport. 2023;(2):4—14. (In Russ.). 

4. Колесников B.M., Мигаль Ю.Ф., Колесников M.B., Сычев A.II., Воропаев А.И. Повышение износостойкости 
тяжелонагруженных трибосистем путем формирования структуры и свойств их контактных поверхностей. Наука 
Юга России. 2022;18(4):59—65. https://www.doi.org/10.7868/825000640220407 

Kolesnikov VI, Migal YuF, Kolesnikov IV, Sitrev AP, Voropaev AP. Increasing the Wear Resistance of Heavily 
Loaded Tribosystems by Forming the Structure and Properties of Their Contact Surfaces. Nauka Yuga Rossii (Science in 
the South Russia). 2022;18(4):59—65. https://www.doi.org/10.7868/S25000640220407 (In Russ.). 

5. Елисеев C.B., Елисеев A.B., Большаков P.C., Хоменко А.П. Методология системного анализа в задачах оценки, 
формирования и управления динамическим состоянием технологических итранспортныхмашин. Москва: «Наука»»; 2021. 679 с. 

Eliseev SV, Eliseev AV, Bolshakov RS, Khomenko AP. Methodology of System Analysis т Problems of Assessment, Formation 
and Control of the Dynamic State of Technological and Transport Machines. Moscow: “Nauka”; 2021. 679 p. (In Russ.). 

6. Ромен Ю.С, Белгородцева Т.М, Дедяев М.В. Состояние ходовых частей вагона и силы взаимодействия в 


системе «экипаж = путь». Транспорт Российской Федерации. 2021;95(4):36-40. | URL: 
https://rostransport.elpub.ru/jour/article/view/127/127 (дата обращения: 19.02.2024). 

Romen YuS, Belgorodtseva TM, Dediaev MV. Condition of Wagon Wheels and Axles and Interaction Forces in the 
"Crew — Track" System. Transport of Russian Federation. 2021;95(4):36—40. https://rostransport.elpub.ru/ 
jour/article/view/127/127 (accessed: 19.02.2024). 

7. Савоськин A.H., Ромен IO.C., Акашев М.Г. Определение вероятностных характеристик боковых сил между 


колесом и рельсом как полезного случайного процесса на фоне помех. Вестник машиностроения. 2022;(4):14—19. 
URL: https://www.mashin.ru/eshop/journals/vestnik_mashinostroeniya/2031/19/ (дата обращения: 19.02.2024). 
Savoskin AN, Romen YuS, Akashev MG. А Useful Random Process of Acting Lateral Forces between а Wheel 
and Rail and [Its Probabilistic Characteristics. Vestnik Mashinostroeniya. 2022;(4):14-19. | URL: 
https://www.mashin.ru/eshop/journals/vestnik_mashinostroeniya/2031/19/ (accessed: 19.02.2024). 
8. Ермоленко И.Ю., Морозов Д.В., Асташков Н.П. Влияние продольных нагрузок на безопасность движения 


при эксплуатации на горно-перевальных участках пути. Вестник Ростовского государственного университета 
путей сообщения. 2021;82(2):104—111. https://doi.org/10.46973/0201-727X_2021 2 104 

Ermolenko IYu, Morozov DV, Astashkov NP. Influence of Longitudinal Loads on Traffic Safety When Operating on 
Mountain Passway Sections. Vestnik RGUPS, Rostov State Transport University. 2021;82(2):104—111. 
https://doi.org/10.46973/0201-727X 2021 2 104 

9. Булдаев A.C., Хишектуева И.Х.Д., Анахин В.Д., Дамбаев Ж.Г. Об одном методе решения задачи 
идентификации динамических систем. Вестник Бурятского государственного университета. Математика, 
информатика. 2020;(4):14—25. https://doi.org/10.18101/2304-5728-2020-4-14-25 


Advanced Engineering Research (Rostov-on-Don). 2024;24(2):125—134. eISSN 2687—1653 


Buldaev AS, Khishektueva IKhD, Anakhin VD, Dambaev ZhG. On One Method for Solving the Problem of 
Identifying Dynamic Systems. Bulletin of Buryat State University. Mathematics, Informatics. 2020;(4):14—25. 
https://doi.org/10.18101/2304-5728-2020-4-14-25 

10. Булдаев A.C. Проекционные методы возмущений в задачах оптимизации управляемых систем. Известия 


Иркутского государственного университета. Серия: Математика. 2014;8:29—43. URL: 
https://mathizv.isu.ru/ru/journal?id=7 (дата обращения: 19.02.2024). 

Buldaev AS. Projection Perturbation Methods in Optimization Problems of Controlled Systems. Bulletin of Irkutsk 
State University. Series: Mathematics. 2014;8:29-43. URL:  https://mathizv.isu.ru/en/journal?id=7 (accessed: 
19.02.2024). 

11. Мижидон А.Д., Хамханов А.К. Гибридная система дифференциальных уравнений, описывающая твердое 


тело, прикрепленное к двум упругим стержням. Вестник Бурятского государственного университета. 
Математика, информатика. 2022;(4):38—47. https://doi.org/10.18101/2304-5728-2022-4-38-47 

Mizhidon AD, Khamkhanov АК. A Hybrid System of Differential Equations Describing a Rigid Body Attached to Two Elastic 
Rods. Bulletin of Buryat State University. Mathematics, Informatics. 2022;(4):38-47. https://doi.org/10.18101/2304-5728- 
2022-4-38-47 

12. Елисеев A.B., Миронов А.С., Елисеев С.В. Формирование математических моделей вибрационных 


взаимодействий элементов технических средств транспортного и технологического назначения. Вестник 
Воронежского государственного университета. Серия: Системный анализ и информационные технологии. 
2022;(1):32—42. https://doi.org/10.17308/sait.2022.1/9199 

Eliseev AV, Mironov AS, Eliseev SV. Formation of Mathematical Models of Vibration Interactions of Elements of 


Technical Means of Transport and Technological Purposes. Proceedings of Voronezh State University. Series: Systems 
Analysis and Information Technologies. 2022;(1):32—42. https://doi.org/10.17308/sait.2022.1/9199 

13. Елисеев A.B., Кузнецов H.K., Елисеев С.В. Новые подходы в оценке динамических свойств колебательных 
структур: частотные функции и связность движений. Труды MAH. 2021;(120):08. https://doi.org/10.34759/trd-202 1- 120-08 

Eliseev AV, Kuznetsov NK, Eliseev SV. New Approaches to the Estimation of Dynamic Properties of Vibrational 
Structures: Frequency Functions and Connectivity of Movements. Trudy МАГ. 2021;(120):08. https://doi.org/10.34759/trd- 
2021-120-08 

14. Елисеев A.B., Кузнецов Н.К. Концепция обобщенного рычага в оценке динамических состояний 


механических колебательных систем в условиях связных вибрационных нагружений. Системы. Методы. 
Технологии. 2023;59(3):7—12. https://doi.org/10.18324/2077-5415-2023-3-7-12 

Eliseev AV, Kuznetsov NK. The Concept of a Generalized Lever in Assessing the Dynamic States of Mechanical 
Oscillatory Systems under Conditions of Connected Vibration Loads. Systems. Methods. Technologies. 2023;59(3):7-12. 
https://doi.org/10.18324/2077-5415-2023-3-7-12 

15. Кашуба B.b., Большаков P.C., Мозалевская A.K., Нгуен Д.Х. Определение реакций связей между 


элементами виброзащитных систем на основе метода структурных преобразований. В: Материалы ХУ 
Всероссийской научно-технической конференции с международным участием «Механики XXI веку». Братск: 
Братский государственный университет; 2016. С. 295—300. 

Kashuba УВ, Bolshakov RS, Mozalevskaya AK, Nguyen Huynh Duc. Identification of Ties Responses between 
Vibroprotection Systems Elements on Base of Structural Transformation Method. In: Proc. XV All-Russian Sci.-Tech. 
Conference with International Participation “Mechanical Engineers to XXI Century”. Bratsk: BrSU; 2016. Р. 295-300. 


Об авторах: 

Роман Сергеевич Большаков, кандидат технических наук, доцент кафедры управления эксплуатационной 
работой Иркутского государственного университета путей сообщения (664074, Российская Федерация, 
г. Иркутск, ул. Чернышевского, 15), ЗРПМ-код: 2025—4049, ORCID, ScopusID, bolshakov_rs@mail.ru 


Валерий Ерофеевич Гозбенко, доктор технических наук, профессор кафедры математики Иркутского 
государственного университета путей сообщения (664074, Российская Федерация, г. Иркутск, 
ул. Чернышевского, 15), 5РІМ-код: 4307—8922, ORCID, ScopusID, vgozbenko@yandex.ru 


Выонг Куанг Чык, соискатель кафедры физики, механики и приборостроения Иркутского государственного 
университета путей сообщения (664074, Российская Федерация, г.Иркутск, ул. Чернышевского, 15), 
ЭРП\-код: 6214—3569, ORCID, ScopusID, trucvq1990@ gmail.com 


Механика 


133 


https://vestnik-donstu.ru 


134 


Большаков Р.С. и др. Возможности оценки динамических состояний транспортных средств 


About the Authors: 
Roman 5. Bolshakov, Cand.Sci. (Eng.), Associate Professor of the Operational Work Management Department, Irkutsk State 


Transport University (15, Chernyshevskogo Str., Irkutsk, 664074, Russian Federation), SPIN-code: 6214—3569, ORCID, 
ScopusID, bolshakov_rs@mail.ru 


Valery E. Gozbenko, Dr.Sci. (Eng.), Professor of the Mathematics Department, Irkutsk State Transport University 
(15, Chernyshevskogo Str., Irkutsk, 664074, Russian Federation), SPIN-code: 4307—8922, ORCID, ScopusID, 


vgozbenko@yandex.ru 


Vuong Quang Truc, postgraduate of the Department of Physics, Mechanics and Instrumentation, Irkutsk State 
Transport University (15, Chernyshevskogo St., Irkutsk, 664074, Russian Federation), SPIN-code: 6214—3569, ORCID, 
ScopusID, trucvg1990@gmail.com 


Заявленный вклад авторов: 

Р.С. Большаков: формирование цели исследования, построение системы уравнений и структурных схем, 
анализ результатов исследования, формулирование заключения. 

В.Е. Гозбенко: корректировка цели исследования, доработка текста, корректировка заключения. 

Выонг Куанг Чык: проведение расчетов, создание графических изображений, анализ и дополнение текста. 


Claimed Contributorship: 
RS Bolshakov: research objective formulation, building a system of equations and block diagrams, analysis of the 


research results, formulation of the conclusion. 
EV Gozbenko: correction of the research objective, revision of the text, correction of the conclusion. 
Vuong Quang Truc: calculation analysis, creation of graphic images, text analysis and addition. 


Конфликт интересов: авторы заявляют об отсутствии конфликта интересов. 
Conflict of Interest Statement: the authors declare no conflict of interest. 

Все авторы прочитали и одобрили окончательный вариант рукописи. 

All authors have read and approved the final manuscript. 


Поступила в редакцию / Received 15.03.2024 
Поступила после рецензирования / Revised 10.04.2024 
Принята к публикации / Accepted 16.04.2024 


Advanced Engineering Research (Rostov-on-Don). 2024;24(2):135—147. eISSN 2687—1653 


МЕХАНИКА 
MECHANICS 


® Check for updates 
BY 


UDC 622.692.4 Original Theoretical Research 
https://doi.org/10.23947/2687-1653-2024-24-2-135-147 


Analysis of the Drag-Reduction Ability of the Layout and Cross-Sectional [=] [=] 

Shapes of Subsea Structures in the Critical Flow Mode ч . 

Henry Е. Annapeh®, Victoria A. Kurushina® 

Industrial University of Tyumen, Tyumen, Russian Federation [=] 
v.kurushina@outlook.com EDN: TDJYXD 

Abstract 


Introduction. Slender structures of subsea energy production systems are under constant influence of currents and waves. 
Hydrodynamic loads result from the interaction of subsea pipelines, umbilicals, equipment supports with fluid flows, and 
lead to the vortex formation in the area behind the structures. Vortex-induced forces are the sources of the cyclic loading. 
They accelerate gradually the fatigue damage, which may result in a failure. One of the ways to reduce the loads on subsea 
structures is to alter the shape of a cross-section, taking into account the flow regime. Dependence of the resulting 
hydrodynamic loads on the cross-sectional shape and relative position of structures has not been studied in details for the 
uniform flow in the critical mode. The current work is aimed at filling this gap. The research objective is to consider the 
impact of the distance between the structures, and also, the presence of a D-shaped structure, placed upstream relative to 
the group of three cylinders of different cross-sectional shapes. 

Materials and Methods. The computational fluid dynamics approach was used in this work for numerical simulations of 
vortex-induced forces in the ANSYS Fluent software for cylinder with D = 0.3 m. Modelling was conducted with the 
Detached Eddy Simulation (DES) method, which combined advantages of the Reynolds-averaged Navier-Stokes equation 
(RANS) method and the Large Eddy Simulation (LES) method. The object of the research was the system of four 
structures in the 2D computational domain, which included the upstream D-shaped cylinder and the main group of three 
cylinders with the circular, squared and diamond shapes of the cross-section. The transient process was considered, where 
structures were under the influence of the uniform flow іп the critical regime at Re = 2.5x10°. 

Results. Five sets of data were obtained in simulations for the time-dependent coefficients of the lift and drag forces: for 
the main system — of the D-shaped, circular, square and diamond structures, and also for the four systems — of only 
D-shaped, only circular, only square and only diamond shaped structures. Additional analysis was conducted for the effect 
of the distance between the structures on the amplitude of fluctuating hydrodynamic force coefficients. The obtained 
results are presented as time histories of coefficients of the lift and drag forces, frequency analysis and contours of 
velocity, pressure and vorticity fields. The results indicate a positive effect of the upstream D-shaped structure on reducing 
the drag force, acting on the central structure in the group of three cylinders located downstream. 

Discussion and Conclusion. The results of the performed studies facilitate the informed decisions regarding the 
arrangement of subsea structures in a group of four objects, depending on the cross-sectional shape and the distance 
between the structures. The upstream D-shaped structure provides reducing the hydrodynamic drag force acting on the 
central structure in the downstream group of three structures, thereby slowing the fatigue accumulation and increasing 
the time of safe operation. 


Keywords: vortex-induced forces, drag coefficient, lift coefficient, uniform flow 
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Оригинальное теоретическое исследование 
Анализ возможности снижения лобового сопротивления за счёт расположения 
и поперечных сечений подводных конструкций в потоке критического режима 


Г.Ф. Аннапе ©, В.А. Kypyumna (54 
v.kurushina(g)outlook.com 


Аннотация 

Введение. Длинные и узкие в поперечнике конструкции морских энергодобывающих систем находятся под по- 
стоянным воздействием течений и волн. Гидродинамические нагрузки являются результатом взаимодействия 
подводных трубопроводов, шлангокабелей, опор оборудования с потоком жидкости и приводят к образованию 
вихрей в зоне за конструкциями. Вихреобразовательные силы служат источником циклического нагружения и 
постепенно ускоряют усталостное разрушение, что может привести к авариям. Одним из способов снижения 
нагрузок на подводные конструкции является изменение формы их поперечного сечения с учетом режима потока. 
Недостаточно изучено, каким образом итоговые гидродинамические нагрузки зависят от формы поперечного се- 
чения и взаимного расположения названных выше элементов систем, находящихся в равномерном критическом 
потоке. Представленная научная работа призвана восполнить этот пробел. Цель исследования — рассмотреть в 
данном контексте значение расстояния между конструкциями, а также наличие полукруглой О-образной кон- 
струкции, размещённой перед группой из трёх цилиндров с разными поперечными сечениями. 

Материалы и методы. Для численного моделирования вихреобразовательных сил использовался метод вычис- 
лительной динамики флюидов в программе ANSYS Fluent для цилиндров диаметром D = 0,3 м. Моделирование 
выполнено методом неприсоединённых вихрей DES, который сочетает в себе преимущества метода усреднён- 
ного по Рейнольдсу уравнения Навье-Стокса RANS и метода крупных вихрей LES. В качестве объекта исследо- 
вания рассматривалась система, состоящая из четырёх конструкций в вычислительном домене в 2D, включая 
стоящий выше по течению полукруглый цилиндр и основную группу из трёх цилиндров круглой, квадратной и 
ромбовидной формы поперечного сечения. Эти конструкции в условиях неустановившегося процесса находятся 
под действием равномерного потока критического режима при Ве = 2,5х10°. 

Результаты исследования. В результате моделирования получены пять наборов данных для изменяющихся во 
времени коэффициентов вихреобразовательных подъёмной силы и силы сопротивления: для основной системы 
из полукруглой, круглой, квадратной и ромбовидной конструкции, а также для четырёх систем из только полу- 
круглых, только круглых, только квадратных и только ромбовидных конструкций. Дополнительно проведён ана- 
лиз влияния расстояния между конструкциями на амплитуду колебаний коэффициентов гидродинамических сил. 
Полученные результаты представлены в виде коэффициентов подъёмной силы и силы сопротивления в дина- 
мике, анализа частот и контуров полей скорости, давления, завихрённости. Результаты позволяют установить 
положительное влияние стоящей выше по течению полукруглой конструкции на снижение силы сопротивления 
на центральную конструкцию в группе из трёх цилиндров ниже по течению. 

Обсуждение и заключение. Результаты проведённых исследований позволяют принимать обоснованные реше- 
ния для расстановки морских конструкций в группе из четырёх объектов в зависимости от формы поперечного 
сечения и расстояния между ними. Установка полукруглой конструкции выше по течению позволяет снизить 
гидродинамическую силу сопротивления на центральную конструкцию в группе из трёх конструкций ниже по 
течению, что замедляет её усталостное разрушение и увеличивает срок эксплуатации. 


Ключевые слова: вихреобразовательные СИЛЫ, коэффициент лобового сопротивления, коэффициент подъемной 
силы, равномерный поток 
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Introduction. Operation and construction of modern offshore systems, specializing on the energy production, 
extraction of resources or carbon capture and storage, require evaluation of the impact of environmental flows on 
equipment and structures. An increased fatigue failure in subsea structures, such as pipelines, risers, cables, piles, 
equipment supports, may come from the vortex shedding phenomenon. The problem is particularly important when 
slender structures are designed to reach deep waters to connect subsystems together. The layout of subsea systems comes 
with the arrangement of structures with different geometry, hydrodynamic properties, and their position in proximity to 
each other. The interference of wakes from these structures and vortex formation patterns is sometimes challenging to 
predict at very high Reynolds numbers due to the turbulent nature of the flow. 

Differences in the flow over a standalone cylinder and two cylinders in tandem are discussed in [1], and three vortex 
shedding regimes are identified for tandem structures. These vortex formation modes include the extended-body regime 
at 1.0 < L/D < 1.8, where L/D is spacing ratio, commonly used to quantify the distance between centres of neighbour 
structures. So that, L corresponds to the distance, and D is the diameter of the structure. In [1], increase of the spacing 
ratio to 1.8 < L/D < 3.8 leads to the reattachment regime, where shear layers detach from the upstream structure and 
reattach to the front side of the downstream structure, so that vortices are formed behind this downstream object. Further 
growth of the spacing ratio, above L/D > 3.8, introduces the co-shedding regime, where a separate vortex is formed from 
the upstream structure and from the downstream structure. Another fundamental research investigated vortex dynamics 
in details through experimental research [2]. One of the following fundamental studies [3] experimentally investigated 
the vortex shedding frequencies of two staggered identical circular cylinders with the Reynolds number Re varying from 
3.2х10* to 7.4x10^, and two fixed side-by-side cylinders at the Reynolds number of 2.5x10* were earlier considered in [4]. 
These investigations provide an important foundation for modern studies in terms of the known effects in fluid forcing 
and vortex shedding patterns. 

Significantly more recent investigations are performed numerically [5] to study the effect of spacing on loads and 
vibrations for two tandem cylinders at subcritical Reynolds numbers, and for specific cases, like a group of mixed large 
and smaller structures [6]. The latter work [6] numerically investigates fluid force coefficients and observes the vortex 
formation pattern on three identical rigid circular cylinders in proximity to a square cylinder. А parametric study is 
conducted in [7] for three identical stationary circular and D-shaped cylinders placed close to a square cylinder at 
Reynolds number 3900 in both linearly and parabolic sheared flows. 

Considering the impact of cross-sectional shapes further, a numerical study is conducted by [8] for a flow over six 
identical stationary cylinders having different cross-sectional shapes at Reynolds number of 2.5x10? in the uniform and 
linearly sheared flow. Rectangular cylinders are investigated in details in [9, 10, 11], where one of the most impactful 
factors for hydrodynamic loads is the aspect ratio of rectangle sides. The works [9, 10] provide new experimental data 
and attempt to develop semi-empirical methods of predicting the response of structures. Further steps in improving the 
modelling approaches for the structural vibration of rectangle-shaped objects under the hydrodynamic excitation are 
performed in [11]. Another branch of studies considers a flow over a sphere [12, 13, 14], while still leaving issues of the 
impact of the cross-sectional shape of subsea structures open. 

Further research on diverse cross-sectional shapes is performed in [15, 16, 17], where triangular and diamond cross- 
sections are studied in comparison. Research [15] is focused on the sensitivity to the corner sharpness for the diamond 
(rhomb) cross-section. Work [16] investigates effects from diverse cross-sections for the system with a rotational degree 
of freedom, when subjected to flow-induced vibration, and study [17] uses cross-sectional shapes for the energy 
harvesting with fluid-structure interaction. 

Following published results, the current work aims to investigate the drag reduction when three structures of different 
cross-sectional shapes are located around a circular cylinder to observe the wake interference and the vortex formation 
pattern between these structures for the spacing ratio L/D varying from 1.67 to 2.83 in the uniform flow at the Reynolds 
number of 2.5x105 with the computational fluid dynamics approach. Specifically, the drag reducing ability of the upstream 
structure is of the research interest. The considered layout is a combination of cylinders placed in tandem, side-by-side 
and staggered position, with a mixture of cross-sectional shapes. 

Materials and Methods. A system of four cylinders with a diameter (characteristic size) of 0.3 m with different cross- 
sectional types is considered in this study. The selected cross-sectional shapes and positions of structures are shown in 
Figure 1. Fluid forces and the flow interference are studied for the spacing ratio of L/D varying from 1.67 to 2.83 in the 
uniform current at the Reynolds number of 2.5x10? and the corresponding velocity of 0.837 m/s. The study focuses on 
the specific layout, where structures are positioned relative to each other in a mixed tandem (cylinders 1 and 3), side-by- 
side (cylinders 2 and 3, and cylinders 3 and 4) and staggered (cylinders 1 and 2, and cylinders 1 and 4) configuration. 
Cross-sectional shapes considered include half-circle, square, circle, and diamond. 
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Periodic 


Shadow 


Fig. 1. Considered structures in the computational fluid domain 


Computational fluid dynamic simulations are performed in the domain with a size of 47D x 27D. The top and bottom 
boundaries are located at a distance of 13.5D away from the center of the circular cylinder, periodic and shadow properties 
are assigned to these boundaries. The left boundary serves as a velocity-inlet, which is located at the distance of 20D from 
the centre of the circular structure 3 in the domain. The value of gauge pressure is set to zero at the pressure-outlet set at 
the right boundary. No-slip conditions are applied to the cylinders. 

The flow around cylinders is simulated using the computational fluid dynamics software ANSYS Fluent, where the 
finite volume method is implemented to solve the Navier-Stokes system. The incompressible flow is considered, and the 
2D DES transient simulations are conducted with the до SST turbulence model. Time integration is performed using 
the second-order implicit transient formulation with a time step of 0.01 s, and the PISO algorithm is used as the solver. 

The DES approach connects capabilities of the Reynolds-Averaged Navier-Stokes (RANS) and Large Eddy 


Simulation (LES) methods [8]. The RANS governing equations for the incompressible flow are as follows: 
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where р is mean pressure, м; is average Cartesian components of the velocity vector, риш; are Reynolds stresses, р is 


density of the fluid and т, is mean viscous stress vector components, which could be expressed as: 


where u is dynamic viscosity. 


(3) 


The Large Eddy Simulation (LES) system of equations for the incompressible flow can be written in the following way: 


where и; and p represent the resolved filtered velocity and pressure, respectively. 
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The diffusion term of the DES model is given by 
Y, = pp'koFprs, 


where В“ is a constant, k stands for fluctuation of the turbulent kinetic energy, œ is specific energy dissipation rate, and 


Еркѕ is as follows: 
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where Caes is a constant, Ajax is local maximum grid map A = (Ai, A», Аз)!. Further, Г, is turbulent length scale: 


JE 


= LCS 8 
F7 pta (8) 
The DES-SST model uses the following zonal formulation: 
L 
Fors -nas| к (1 Бы). (9) 
des тах 


where Fssr = 0, Ру, Р, and Ру, F are mixed functions of the SST model. 

Table 1 below provides results of the mesh independence test for the uniform flow of Re = 2.5х10°. The mesh settings 
are adopted from [7], and the accuracy of the grid is demonstrated by comparisons in Table 1. All subsequent analysis in 
this paper is performed with Mesh 2, and the results include signals and frequencies of the fluid force coefficients and an 
indication of the vortex shedding pattern. The drag force fluctuations are presented in this work in terms of the drag force 
coefficient Cp which comprises the mean drag coefficient Cpo and the fluctuating drag coefficient Cp", as follows: 


Cp =C + Ch. (10) 
The lift force fluctuations are presented using the lift coefficient Cz. 
Table 1 
Mesh independence test results 
Cases Cpo Number of cells Strouhal number 
Ве = 2.5x10° 
Current study 
Mesh 1 0.98 63,345 0.24 
Mesh 2 1.08 86,478 0.24 
Mesh 3 1.08 131,041 0.24 
Published data 
Lehmkuhl, et al. (2014) (LES) [18] 0.833 = 0.238 
Achenbach&Heinecke (1981) (Experiment) [19] 1.135 = 0.230 
Re = 3,900 
Current study, Mesh 2 0.93 86,478 0.18 
Wornom, et al. (2011) (VMS-LES) [20] 0.99 E Е 
Ве = 3.6x106 
Current study, Mesh 2 0.4100 86,478 - 
Porteous, et al. (2015) (URANS) [21] 0.4206 = – 
Nazvanova, et al. (2022) (URANS) [22] 0.4657 74,496 = 


Results. In this study, simulations are performed in two series. The first series is focused on recognising the overall 
effect of various cross-sectional shapes, placed at L/D = 2.00 from each other. Cylinder numbers here correspond to the 
ones used in Figure 1, according to the cylinders position. The calculation results for this set are reported in Table 2, in 
comparison to the case of structures with mixed cross-sections. The second series provides an insight into the impact of 
L/D ratio on hydrodynamic loads observed for the mixed cross-section case only, as in Figure 1. These results are 
summarised in Table 3 and allow defining the drag reducing effect of the D-shape upstream structure on loads when 
cross-sections are different. 
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Simulation results for the same arrangement with different cross-sectional shapes for L/D = 2.00 


Table 2 


LID =2.0 Em All circular | АП square All D-shaped All diamond-shaped 
Giossasctionkcin Fig. d structures structures structures structures 
Cylinder 1 
Cpo 0.45 0.28 0.49 0.48 0.46 
cf 0.19 0.14 0.43 0.28 0.17 
CL 0.18 0.04 0.07 0.13 0.28 
Cylinder 2 
Сро 0.87 0.39 0.89 0.64 0.70 
СЇ 1.03 0.38 1.38 0.49 0.60 
CL 1.38 0.53 0.07 0.11 0.07 
Cylinder 3 
Coo 0.3 0.24 0.26 0.35 0.70 
cf 0.45 0.37 0.91 0.68 0.63 
CL 1.03 0.73 1.26 0.30 0.73 
Cylinder 4 
Сро 0.98 0.39 0.83 0.61 0.99 
СЇ 0.90 0.27 0.96 0.44 1.05 
Ст 1.22 0.81 1.63 0.37 1.11 
Table 3 
Simulation results for the mixed cross-sections at various L/D ratio 
Cylinder 1 Cylinder 2 Cylinder 3 Cylinder 4 
L/D 
C» | Ch Ci C» | Ch Ci C» | Ch Ci C» | Ch Ci 
1.67 0.14 0.18 0.15 0.98 1.17 1.44 0.48 0.57 0.66 0.93 0.75 1.13 
1.83 0.33 0.21 0.18 0.93 1.54 1.40 0.33 0.50 0.86 1.01 1.36 1.15 
2.00 0.45 0.19 0.18 0.87 1.03 1.38 0.30 0.45 1.03 0.98 0.90 1.22 
2.17 0.49 0.05 0.19 0.81 0.14 1.11 0.25 0.06 0.83 0.96 0.38 1.10 
2.33 0.37 0.12 0.00 0.78 0.34 1.03 0.25 0.06 0.28 1.00 0.07 0.79 
2.50 0.57 0.16 0.20 0.76 0.89 1.18 0.26 0.47 0.77 1.01 0.64 1.06 
2.67 0.57 0.02 0.27 0.78 0.16 1.37 0.28 0.08 0.82 1.02 0.03 0.96 
2.83 0.58 0.21 0.21 0.84 0.91 1.70 0.20 0.52 1.16 0.93 0.68 1.03 


Comparison of all circular, all square, mixed (as in Fig. 1), all diamond and all D-shapes with each other in Table 2 
reveals relatively lower hydrodynamic loads for all four structures observed for the circular shapes at the L/D = 2.00. The 
largest mean drag coefficient here is experienced by the fourth structure in both mixed-shaped and all diamond-shaped 
arrangement. The highest maximum fluctuating drag coefficient of 1.38 is observed for cylinder 2 in the all square-shaped 
arrangement. The largest maximum amplitude of the lift coefficient of 1.63 also belongs to the all square-shaped 
arrangement, but corresponds to cylinder 4. Mixed-shaped arrangement (or basic case of structures with alternate cross- 
sections, as in Fig. 1) at L/D — 2.00 demonstrates a relative consistency in large amplitudes of the lift coefficient for 
cylinders 2, 3, 4, which makes estimations of hydrodynamic loads for this arrangement more important, due to higher 


expected loads. 


л 
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Fig. 2. Time histories of fluctuating drag coefficients for four cylinders of a different cross-sectional shape with L/D — 2.0: 
a — cylinder 1; b — cylinder 2; c — cylinder 3; d — cylinder 4 
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Fig. 3. Time histories of lift coefficients for four cylinders of a different cross-sectional shape with L/D — 2.0: 
a — cylinder 1; b — cylinder 2; c — cylinder 3; d — cylinder 4 


Further observation of signals of the fluctuating drag and lift coefficient in Figures 2 and 3 reveals meaningful instabilities 
in forces experienced by all structures with square shapes, and for structures 2, 3, 4 with diamond and mixed shapes. The 
comparison shows that circular and D-shaped structures would experience lower and more stable fluid loads in this arrangement. 
Table 2 and Figures 2 and 3 confirm that cylinder 1 in the shape of a circle or D-shape provides reducing fluid loads for the 
three downstream cylinders. This substantiates the common interest in further exploration ofthe effect of the D-shaped structure 
on reducing the drag force in the arrangement with all alternative cross-sectional shapes. 
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This effect is studied in details in the next (second) simulation series, presented for each cylinder in Figures 4—7 in 
terms of the fluid forces and in Figures 8, 9 — in terms of the fluid flow characteristics for the considered computational 
domain. Figure 4 illustrates fluid loads on the upstream structure, where the most unstable signal (at L/D = 2.17) belongs 
to fluctuations of the drag force. Figure 4 c also indicates presence of multiple frequencies in signals of the fluctuating 
drag coefficient, while a single dominating frequency can be identified for the lift force coefficient in Figure 4 d. 
Figures 4 c and 4 d demonstrate that the frequency of both lift and drag forces generally increases with the growing L/D 
ratio, and the maximum frequency is indicated by signals at L/D = 2.67 and 2.83. 
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Fig. 4. Time histories of fluid force coefficients for cylinder 1 in the uniform flow: 
a — fluctuating drag coefficient; b — lift coefficient; c — fluctuating drag coefficient FFT; d — lift coefficient FFT 


Relatively similar complexity of frequencies of the lift force is observed in Figure 6 d for cylinder 3, where one to 
two dominant frequencies could be clearly identified. At the same time, more than two frequencies are observed in 
Figure 6 c for each signal of the fluctuating drag force. The pattern of growth in the overall dominant frequency with the 
increasing L/D is still recognisable for cylinder 3, similar to cylinder 1. Combination of frequencies is even more complex 
for cylinders 2 and 4, as comes from Figures 5 c-d and 7 c-d, this is not possible to indicate clear dependences in the 
frequencies of the fluctuating drag force. Some resemblance of the found growth trend could be still observed in Figure 5 d 
for the lift force coefficient for cylinder 2. This provides evidence for the generally unstable nature of hydrodynamic 
loads acting on the square cross-section shown in Figure 5. 

Impact of the L/D ratio on the mean drag coefficient for cylinder 1, based on Table 3, is partially confirmed: apart 
from a couple of deviations, the mean drag force increases with the growth of L/D. Data for cylinders 2 and 4 do not 
indicate a specific pattern, as the values fluctuate back and forth within 10% from the initial mean drag coefficient at the 
smallest L/D. Cylinder 3, on the contrary, indicates a stable pattern of the reduced mean drag coefficient with the increased 
L/D, so that the reducing ability of the upstream D-shaped structure is evident, but for the central cylinder only. The 
largest mean drag coefficient of 1.02 is linked to cylinder 4 at L/D = 2.67. The considered range of L/D allows observing 
an important transition from the strong to minor interference in the wake of the three paired structures. 

The highest fluctuating drag coefficients of 1.54 and 1.36 are linked to cylinders 2 and 4, respectively, both observed 
at L/D = 1.83. The feature of the maximum amplitude of the fluctuating drag coefficient, indicated in Table 2, is in absence 
of a specific dependence from the L/D, the values rapidly change from near zero to relatively high with a small increment 
of change in the ratio. The fluctuating drag coefficient has generally the lowest amplitudes for cylinder 1, average 
amplitudes — for cylinder 3, and the largest amplitudes — for cylinders 2 and 4. 

The largest maximum amplitude of the lift coefficient occurs at L/D = 2.83 for cylinder 2. The lift coefficient 
generally resembles the distribution, similar to the maximum amplitude found for the fluctuating drag coefficient: the lift 
force appears to be the smallest for cylinder 1, relatively intermediate — for cylinder 3, and the highest — for cylinders 
2 and 4, with no specific pattern linked to the L/D increase and reduction. This allows us to conclude that the ability to 
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reduce hydrodynamic forces by placing the upstream D-shaped structure in front of the array is limited. The force 
reduction is observed mainly for the structure placed in tandem downstream, and the effect is most pronounced for the 
mean drag coefficient, with some reduced effects also seen for the fluctuating drag and lift coefficients. 
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Fig. 5. Time histories of fluid force coefficients for cylinder 2 in the uniform flow: 
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Fig. 6. Time histories of fluid force coefficients for cylinder 3 in the uniform flow: 
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Fig. 7. Time histories of fluid force coefficients for cylinder 4 in the uniform flow: 


c) 


Vorticity 
70.0 
56.0 
42.0 
28.0 
14.3 


0.0 


d) 


Fig. 8. Contours of the flow characteristics for L/D — 1.67 at 200 seconds: 
a — velocity contour; b — vorticity contour; c — pressure contour; d — streamline 
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с) а) 


Fig. 9. Contours of the flow characteristics for L/D = 2.83 at 200 seconds: 
a — velocity contour; b — vorticity contour; c — pressure contour; d — streamline 


Figures 8—9 show the velocity, vorticity, pressure and streamlines of the flow around cylinders for some selected L/D, 
where both proximity and wake interference among the cylinders are presented for the time step of 200 s. The flow around 
cylinders is complex, and vortex formation patterns are highly affected, as the distance between the cylinders increases. 
The proximity interference is observed for cylinders 2, 3, and 4, alternate single vortices are shed on the downstream side 
of these structures. For the wake interference at L/D = 1.67, free shear layers separate from the upstream cylinder 1 and 
reattach themselves to the upstream side of cylinder 3, and a vortex street is only formed at the downstream side of 
cylinder 3. At this distance, a broad region of wake is created at the downstream side of cylinders 2, 3, and 4. As the L/D 
increases above 2.00, there are vortices formed at the upstream side of cylinder 3, in the wake of cylinder 1. A vortex 
street is also formed at the downstream side of cylinder 3, with formation of 2S vortices. Figures 8 a, 9 c demonstrate a 
group of minor vortices formed following the diamond-shaped cylinder 4 and a vortex pair formed in a similar to S+P 
vortex cycle past cylinder 2. 

Discussion and Conclusion. The 2D numerical simulations are performed in this work for cylinders with different 
cross-sectional shapes at the Reynolds number of 2.5x10? using the DES approach. The considered cylinders are studied 
in a complex position of an upstream D-shaped structure in front of three paired structures, with the aim to investigate the 
drag-reducing ability of this specific layout, observe the flow complexity, the wake interference from each structure, and 
vortex formation patterns. 

The following conclusions could be made from this study: 

1. The ability to reduce hydrodynamic forces by placing the upstream D-shaped structure 15 mainly limited to the 
structure placed in tandem downstream, and the effect is most pronounced for the mean drag coefficient. 

2. Overall, the mean drag coefficient of cylinders is observed to be affected by varying L/D, with the main effect on the 
mean drag coefficient of cylinder 1, which grows with increasing L/D, and of cylinder 3, which reduces with increasing L/D. 

3. Competition of frequencies is observed for the fluctuating drag coefficient for all structures and lift coefficient 
signal of cylinders 2 and 4. This competition is due to the joint effects of both the uniform current and wake interference, 
which intensifies at a lower L/D in terms of changes to the resulting vortex street. 

4. Both proximity and wake interference among the cylinders are observed. The flow around cylinders is complex, 
and vortex formation patterns are highly affected as the distance between the cylinders increased with 2S being the major 
vortex type formed and shedding the additional vortices from the square and diamond structures. 
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The study generally contributes to the field of knowledge by advancing our understanding of fluid-structure 
interactions, drag reduction strategies, and vortex dynamics, with potential applications in offshore energy systems. The 
current work contributes to the development of the drag reduction strategies through analyzing the impact of the upstream 
D-shaped structure on downstream cylinders. Understanding how different structural configurations affect drag can 
inform the design of more efficient systems in various engineering applications, such as subsea transportation of fluids. 
By observing flow complexity, wake interference, and vortex formation patterns, this study contributes to the 
understanding of fluid dynamics around complex geometries. This knowledge is crucial for optimizing the performance 
of structures in environments where fluid flow plays a significant role, such as subsea engineering. The results of the 
present research highlight the effect of varying the aspect ratio L/D on drag coefficients to inform engineering designs of 
similar arrangements. This study reveals the intricate vortex dynamics and shedding patterns, particularly concerning the 
proximity and wake interference among cylinders. Understanding these phenomena can aid in predicting and controlling 
flow behavior around complex configurations, leading to more efficient designs and better performance in practical 
applications in offshore systems. 
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Аннотация 

Введение. Устройства сбора и накопления энергии из внешней среды представляют собой маломощные источ- 
ники электрической энергии, которые активно используются, в том числе в автономных приборах мониторинга 
поврежденного состояния различных конструкций. Рабочим элементом этих устройств является пьезоэлектриче- 
ский генератор (ПЭГ) — преобразователь механической энергии в электрическую. Конструирование ПЭГ свя- 
зано с предварительным построением их математических и компьютерных моделей, с помощью которых произ- 
водится расчет и оптимизация конструкций. Одним из способов моделирования и расчета ПЭГ является разра- 
ботка приближенных методов расчета на основе прикладных теорий. В литературе известны и ранее разработаны 
прикладные теории расчета изгибных колебаний многослойных пьезоактивных пластин. Однако информации об 
изгибно-сдвиговых колебаниях, как инструменте повышения эффективности инженерных расчетов описанных 
конструкций, в научной литературе недостаточно. Целью настоящей работы являлась разработка прикладного 
метода расчета изгибных и сдвиговых колебаний пьезокерамических пластин, в том числе пористых. 
Материалы и методы. В качестве пьезоактивного материала пластины используется пьезокерамика PZT-4, в 
том числе пористая. При использовании пористой керамики жесткость конструкции уменьшается в большей сте- 
пени, чем пьезомодули, что позволяет получить более эффективный ПЭГ при механическом воздействии. Мате- 
матическая постановка осуществлена в рамках линейной теории электроупругости при поляризации пластины 
по толщине. Боковые стороны пластины электродированы, правая сторона закреплена, а на левой задан гладкий 
контакт в вертикальной стенке. Установившиеся колебания пластины вызываются давлением на лицевые поверх- 
ности пластины или разностью электрических потенциалов на электродах. Для расчета характеристик ПЭГ в ра- 
боте предлагается прикладная теория, основанная на гипотезах о распределении характеристик напряженно-де- 
формированного состояния и электрического поля. 

Результаты исследования. Рассмотрены поперечные колебания пьезокерамической пластины в низкочастотной 
области (ниже первого изгибно-сдвигового резонанса). В силу того, что математическая постановка рассмотрена 
в рамках линейной теории упругости, задача разделилась на сумму двух. В первой учитывалось механическое 
воздействие: на лицевые поверхности пластины действует распределенная нагрузка и поперечная сила на левом 
конце, а потенциалы на электродах равны нулю. Во второй задаче механические нагрузки отсутствовали, но за- 
давалась разность потенциалов на электродах. На основе гипотез о распределении деформаций, механических 
напряжений и электрического потенциала обе задачи были сведены к системе обыкновенных дифференциальных 
уравнений и граничных условий. Сравнение с результатами расчетов методом конечных элементов в пакете 
АСЕГАМ показали адекватность предложенной прикладной теории в низкочастотной области. 


© Соловьев A.H., Чебаненко B.A., Оганесян П.А., Фоменко E.H., 2024 


Соловьев А.Н. и др. Об одном методе расчета изгибных и сдвиговых колебаний пористого пьезоэлемента в низкочастотной области 


Обсуждение и заключение. Поскольку постановка задачи рассматривалась в линейной теории электроупругости и 
изучалась низкочастотная область, в работе удалось задачу об изгибных и сдвиговых колебаниях пластины из по- 
ристой пьезокерамики разделить на две: изгибную — с механическим воздействием при нулевых потенциалах и 
сдвиговую — при задании разности потенциалов и нулевом механическом воздействии. Использованы соответ- 
ствующие гипотезы об изгибе и сдвиге, построены две системы обыкновенных дифференциальных уравнений и 
граничных условий, которые решаются аналитически без использования «тяжелых» конечно-элементных пакетов. 
Для сравнения результатов и подтверждения адекватности предложенного метода проведено конечно-элементное 
моделирование таких задач B специализированном пакете АСЕГАМ. Это сравнение показало, что ошибка в опреде- 
лении смещений и электрического потенциала при использовании этого подхода, в случае задания механических 
нагрузок и разности потенциалов, не превышает 6 %. Разработанный в статье метод может быть применен при про- 
ектировании пьезоэлектрических генераторов накопления энергии в низкочастотной области. 
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пластины, сдвиг пластины, прикладная теория 
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Abstract 

Introduction. Devices for collecting and storing energy from the external environment are low-power sources of electric 
energy that are actively used. The autonomous devices for monitoring the damaged condition of various structures include 
them as well. The working element of these devices is a piezoelectric generator (PEG) — a converter of mechanical 
energy into electrical energy. The design of PEG is associated with the preliminary construction of their mathematical 
and computer models, with the help of which the calculation and optimization of structures is carried out. One of the ways 
to model and calculate PEG is to develop approximate calculation methods based on applied theories. The applied theories 
for calculating bending vibrations of multilayer piezoactive plates are known and previously developed in the literature. 
However, in the scientific literature there is not enough information about bending and shear vibrations as a tool for 
improving the efficiency of engineering calculations of the described structures. The objective of this work was to develop 
an applied method for calculating bending and shear vibrations of piezoceramic plates, including porous ones. 

Materials and Methods. Piezoceramics PZT-4, including porous ones, were used as the piezoactive material of the plate. 
When using porous ceramics, the rigidity of the structure decreased to a greater extent than the piezoelectric modules, 
which made it possible to obtain a more effective PEG under mechanical action. The mathematical formulation was 
carried out within the framework of the linear theory of electroelasticity with plate polarization in thickness. The sides of 
the plate were electrodated, the right side was fixed, and a smooth contact in the vertical wall was set on the left side. 
Steady-state vibrations of the plate were caused by pressure on the front surfaces of the plate or the difference in electrical 
potentials at the electrodes. To calculate the characteristics of PEG, the authors proposed an applied theory based on 
hypotheses about the distribution of characteristics of the stress-strain state and the electric field. 

Results. Transverse vibrations of a piezoceramic plate in the low-frequency region (below the first bending-shear 
resonance) were studied. Due to the fact that the mathematical formulation was considered within the framework of the 
linear theory of elasticity, the problem was divided into the sum of two. The first one took into account the mechanical 
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effect: a distributed load and a transverse force at the left end acted on the front surfaces of the plate, and the potentials 
at the electrodes were zero. In the second task, there were no mechanical loads, but the potential difference was set at the 
electrodes. Based on hypotheses about the distribution of deformations, mechanical stresses and electric potential, both 
problems were reduced to a system of ordinary differential equations and boundary conditions. Comparison with the 
results of calculations by the finite element method in the ACELAN package showed the adequacy of the proposed applied 
theory in the low-frequency region. 

Discussion and Conclusion. Since the formulation of the problem was considered in the linear theory of electroelasticity, 
and the low-frequency region was studied, the work succeeded in dividing the problem of bending-shear vibrations of a 
porous piezoceramic plate into two: bending — with mechanical action at zero potentials, and shear — when setting the 
potential difference and zero mechanical action. The corresponding hypotheses about bending and shear were used. Two 
systems of ordinary differential equations and boundary conditions, which were solved analytically without the use of 
“heavy” finite element packages, were constructed. To compare the results and confirm the adequacy of the proposed 
method, the finite element modeling of such tasks was carried out in a specialized ACELAN package. The comparison 
showed that the error in determining displacements and electric potential when using this approach, in the case of setting 
mechanical loads and potential differences, did not exceed 6%. The method developed in the paper can be applied in the 
design of piezoelectric generators for energy storage in the low-frequency region. 


Keywords: energy collection device, piezoelectric generator, porous ceramics, plate bending, plate shear, applied theory 
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Введение. Пьезоэлектрические генераторы (ПЭГ) используются для преобразования механической энергии в 
электрическую с последующим ее накоплением. Одна из областей применения ПЭГ — создание маломощных 
автономных возобновляемых источников электрической энергии. Рабочим элементом ПЭГ является 
пьезокерамический элемент определенной формы. Форма и тип деформации этого элемента определяют 
пьезомодуль, который характеризует преобразование механической энергии деформации в электрическую. Так 
пьезомодуль 233 связан с растяжением-сжатием вдоль оси поляризации, з — c такой же деформацией в поперечном 
направлении к этой оси, 415 — со сдвигом. Использование пористой керамики позволяет создавать более 
эффективные ПЭГ. Это связано с тем, что модули упругости пористой керамики с ростом пористости убывают 
значительно сильнее, чем пьезомодули. Таким образом при одной и той же механической нагрузке амплитуда 
деформации у пористой керамики будет больше, следовательно, выходной электрический потенциал также больше. 

Расчет ПЭГ может быть произведен методом конечных элементов, реализованном в пакетах ANSYS, 
ACELAN, COMSOL и других. Для пьезоэлементов, один или два размера которых значительно меньше других 
(пластины, стержни), могут быть построены прикладные теории расчета на основе гипотез о распределении ме- 
ханического и электрического полей. Без применения «тяжелых» конечно-элементных пакетов прикладные тео- 
рии позволяют моделировать различные устройств на основе пьезоактивных материалов. В качестве таких мате- 
риалов рассматривались пьезоэлектрические, пьезомагнитные и композиционные пьезомагнитоэлектрические. 
Построение этих теорий основано на принятии гипотез о распределении механических, электрических и магнит- 
ных полей. Эти гипотезы связаны с характером колебаний упругих и пьезоактивных элементов ПЭГ. Наиболее 
распространенными конструкциями являются активные и полупассивные биморфы на основе многослойных пла- 
CTHH, поляризованные по толщине с электродами на лицевых поверхностях, совершающие поперечные изгибные 
колебания. Исследованию устройств со сдвиговой деформацией пьезоэлементов посвящен ряд работ. Электри- 
ческая модель с пьезоэлектрическими определяющими уравнениями режима dis и модель с одной степенью сво- 
боды объединены для описания характеристик сбора энергии пьезоэлектрического кантилевера в сдвиговом ре- 
жиме в работе [1]. Предлагаемая модель используется для моделирования частотной зависимости выходного пи- 
кового напряжения и мощности. Результаты показывают хорошее совпадение с экспериментом и конечно-эле- 
ментным расчетом в ANSYS. В работе [2] разработан пьезоэлектрический преобразователь энергии сдвигового 
режима для использования энергии потока воды под давлением. Он преобразует энергию потока в электрическую 
энергию путем пьезоэлектрического преобразования с колебанием пьезоэлектрической пленки. Разработана мо- 
дель конечных элементов для оценки генерируемого напряжения пьезоэлектрической пленки, которая хорошо 
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согласуется с проведенным натурным экспериментом. Одномерная полностью связанная модель колебания 
балки на основе гипотез типа Тимошенко, которая предоставляет единую общую основу для анализа энергии в 
режиме сдвига и изгиба, представлена в работе [3]. В работе [4] изучается влияние неоднородности свойств пла- 
стины при сдвиговых и крутильных колебаниях ее центральной части. В экспериментальной работе [5] был пред- 
ставлен многослойно-цилиндрический пьезоэлектрический сдвиговый актуатор (MCPSA), работающий в режиме 
сдвига dis, для прецизионного срабатывания при большой механической нагрузке. Актуатор был изготовлен из 
пьезоэлектрических керамических колец Pb(Zr,Ti)O3 (PZT-51), которые были концентрически собраны вместе B 
электрически параллельном соединении с попеременно положительной и отрицательной поляризацией в осевом 
направлении. В работе [6] создан метаматериал из идентичных элементарных ячеек, спроектирован и изготовлен 
искусственный прототип устройства с характерными узорчатыми электродами и расположенными в ряд пьезо- 
керамическими субъединицами, который, как доказано, идеально генерирует синтетическую деформацию сдвига 
грани. При том же напряжении возбуждения наблюдается усиление смещения сдвигового типа более чем на по- 
рядок, по сравнению с предыдущими объемными элементами в режиме dis. В статической постановке в работе [7] 
теоретически установлено поле электромеханической связи в сдвигово-изгибающем режиме для кольцеобразной 
пьезоэлектрической пластины. В соответствии с классической теорией упругих пластин малого изгиба и пьезо- 
электрическими определяющими уравнениями было достигнуто аналитическое решение изгибной деформации 
пьезоактюатора под действием электрического поля и концентрированной или равномерно распределенной ме- 
ханической нагрузки. Механизм создания изгибной деформации объясняется осесимметричной сдвиговой де- 
формацией, которая дополнительно вызывает изгибную деформацию одной пьезоэлектрической пластины B 
форме кольца. Этот механизм существенно отличается от механизма пьезоэлектрических биморфных или уни- 
морфных приводов, о которых сообщалось ранее. Проведена оптимизация конструкции кольцеобразного пьезо- 
актуатора. В работе [8] на основе одномерной модели строится функция отклика датчика на основе сдвиговых 
резонаторов (срезы кварца) объемной акустической волны, которые перспективны для поточных измерений вяз- 
кости жидкости, например, в промышленных процессах. В работе [9] с помощью метода конечных элементов 
исследовался пьезоэлектрический преобразователь управления высотой полета с использованием деформации 
модели сдвига. В [10] теория функционально-градиентной пластины с четырьмя неизвестными сдвиговой дефор- 
мации применяется для выражения компоненты смещения. Распределение электрического потенциала представ- 
ляет собой линейную функцию по толщине. Пластина находится под механической нагрузкой и электрическим 
напряжением. Основные уравнения и граничные условия выводятся с использованием принципа виртуальной 
работы. Проведен анализ напряжений и деформаций от параметров конструкции. Электромеханический анализ 
потери устойчивости пьезоэлектрической нанопластины при сдвиге с использованием модифицированной тео- 
рии парных напряжений с различными граничными условиями изучался в работе [11]. Чтобы учесть электриче- 
ские эффекты, к пьезоэлектрической нанопластине прикладывали внешнее электрическое напряжение. Была ис- 
пользована упрощенная теория сдвиговой деформации первого порядка. Основные дифференциальные уравне- 
ния были получены с использованием принципа Гамильтона и нелинейных деформаций Фон-Кармана. В итоге 
результаты показали, что влияние внешнего электрического напряжения на критическую сдвиговую нагрузку, 
возникающую на пьезоэлектрической нанопластине, незначительно. В работе [12] с помощью комбинации двух 
классических подходов моделирования нелинейного поведения пьезоэлектрических материалов исследуется пье- 
зоэлектрический привод сдвигового типа для атомно-силового микроскопа. В частности, новизна предлагаемого 
метода состоит в том, что он сочетает в себе два источника нелинейности полезависимой модели Мюллера и 
Чжана с частотно-зависимой моделью Дамьяновича. Численные результаты, полученные с помощью метода ко- 
нечных элементов (МКЭ), сравниваются с экспериментом. 

Менее исследованы в научной литературе колебания, в которых кроме изгиба присутствует сдвиг, т.е. «рабо- 
чим» является пьезомодуль dıs, значение которого убывает с увеличением пористости, HO в меньшей степени, 
чем упругие модули. Последнее обстоятельство позволяет построить эффективное устройство преобразования 
энергии. Поэтому разработка прикладной теории расчета ПЭГ с использованием пористой пьезокерамики на ос- 
нове упрощенных моделей без использования «тяжелых» конечно-элементных пакетов представляется весьма 
актуальной задачей. Целью настоящей работы явилось построение прикладного метода расчета для установив- 
шихся поперечных колебаний в низкочастотной области пористой пьезокерамической пластины, характеризую- 
щиеся как сдвигом, так и изгибом. 

Материалы и методы. Исследуемый ПЭГ представляет собой пьезокерамическую пластину (длина [, 
толщина Å), поляризованную по толщине, консольно-закрепленную по правой боковой стороне, левая боковая 
сторона прикреплена к инерционной массе, которая совершает вертикальные колебания и закреплена в 
горизонтальном направлении. Электроды расположены на боковых сторонах пластины, поэтому при разности 
потенциалов на них и отсутствии механической нагрузки, основной деформацией в низкочастотной области 
является сдвиг. Рассматриваются колебания, частота которых меньше частоты первого резонанса. 
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Математическая постановка задачи 
Математическая постановка задачи описывается системой дифференциальных уравнений [13] и соответству- 
ющими граничными условиями. 
PpUtagpju-V-o=f; V-D=0 
с=с? (+В) е^ -E Е = -УФф (1) 
р+с̧ад =е,-(&+642)+э7-Е e=(Vu+Vu")/2 
При рассмотрении пористой керамики связности 3-0 в уравнениях (1) используются эффективные физиче- 


ские константы, определенные с помощью пакета ACELAN-COMPOS [14]. Эти эффективные свойства на основе 
представительных объемов (рис. 1) получены в работе [15] и представлены в таблице 1. 


Рис. 1. Представительные объемы в пакете ACELAN-COMPOS [14] 


Таблица 1 
Эффективные свойства пористой керамики 

% пористости 0 10 20 30 40 50 60 70 80 
р, kg/m? 7500 | 6750 | 6000 | 5250 | 4500 | 3750 | 3000 | 2250 | 1500 
cE 1010, N/m? 13,90 | 11,56 | 9,25 6,85 5,05 3,34 2,07 1,26 0,68 
cE 1010, N/m? 7,78 6,15 4,66 3,14 | 2,10 1,16 0,62 0,28 0,13 
ch, 1019, N/m? 7,43 5,82 4,25 2,82 1,87 1,06 0,52 0,24 0,10 
che 1010. N/m? 11,50 | 9,53 7,23 5,42 3,91 2,72 1,63 0,91 0,47 
che 1010, N/m? 2,56 2,23 1,83 1,44 1,10 0,74 0,44 0,23 0,10 
es , Chm? 15,10 | 13,38 | 11,37 | 9,59 7,68 5,93 3,93 2,30 1,25 
es, Chm? -5,20 | -4,23 | -3,14 | -2,07 | -1,32 | -0,75 | -0,43 | -0,21 | -0,10 
est , Chm? 12,70 | 10,96 | 8,96 6,91 5,00 3,30 1,95 100 | 0,44 
ir gn 730 663 582 509 439 349 263 191 122 

ке / ео 635 567 492 413 345 270 199 130 75 


По данным таблицы 1 построены зависимости эффективных упругих модулей и пьезомодулей от процента 
пористости [15], которые представлены на рис. 2. 


100 500 
<) & 300 ; 
Е 60 = Ef ef gs 200 = ай 
E — EET Е — 49 
a 40 3 3 100 — ad 
d D 

+ * 
p & 100 Lo 
0 Е 
0 20 40 60 80 "m 30 40 50 60 70 80 
Porosity, 96 Porosity, 96 
a) 6) 


Рис. 2. Зависимости значений or процента пористости [15]: а — упругих моделей; 6 — пьезомодулей 
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В соответствии с результатами, представленными на рис. 2, модуль упругости Е! убывает значительно быст- 
рей, чем пьезомодуль 415. 

Построение прикладной теории 
Построенный метод расчета состоит из решения двух задач: в первой рассматривается изгиб под действием ме- 
ханической нагрузки при нулевой разности потенциалов, во второй — моделируется сдвиг, вызванный разностью 
потенциалов при отсутствии механической нагрузки. В обоих случаях учитывается отсутствие зарядов на лице- 
вых поверхностях пластины. При действии механической нагрузки и разности электрических потенциалов ре- 
зультаты двух задач складываются в силу линейности задачи. 

Действие механической нагрузки, потенциалы на электродах равны нулю. Принимаются гипотезы типа 
Кирхгофа-Лява относительно равенства нулю нормальных напряжений, а для перемещений: 


а 
из =02> (x), ш =-| — 022 (х) |z. (2) 
dx 
Распределение электрического потенциала по толщине предполагается квадратичным и симметричным: 
2 2 22 
ф = Ф| (22-1) + Ф (1+ ie t Ф, (х) Its Б (3) 
Для учета граничных условий на концах пластины (x = 0, /) получено выражение для поперечной силы: 


8сззозз +8ез 8 


C13€33 2 p? 
2 
1 (cgs +e) 8033233 + 8e33 8 h? 
31 


Qı +e 
12 C33 (сззазз +e )h? he 


1 2 
В 1 8c33g33 + 8е33 4 3 disi (Zo) 


12 13 (casgas ei) h А? dx T 
e33 (свеззй? — сзземй?) 
C13 A 3 C13 А А 
1 (casgas +e )h exi(eisessh —c33e31h ) 
+ C11 + 
12 C33 (casg33 +033) h? 
2 2 
1 es (cressh — Сззезій )^ d? 
2 7022 (х) 
24 C33 233 + €33 dx 


С учетом равенства нулю нормальной компоненты вектора электрической индукции Ha лицевых поверхно- 
стях (Z = +h), уравнения для неизвестного прогиба UZ»(x) и распределения электрического потенциала D(x) 
имеют вид: 


8сззезз +8ез 8 


C13€33 


(c +e? ) h’ 2 
1 33833 33 ion 8033233 + 8e33 8 hà 
31 
2 C33 (casgss teh)? h? 
1 | 8c33g33 +8e33 4 | 3 í 
— — 615 €f. Эа 12 h —esh —— 5 (x) + 
2 2 (casgas +е%)л h dx (5) 
€33 (свеззй? сузем?) 
C13 2155 C13 j 4 
1 (c33833 +e )h €31 (cisessh — сззезіћ ) 3 
+ C11 + 
2 C33 (casgas cen? 


] 815 сзеззй” — esses i? h |( q^ 
l ) z UZa (x) -WphUZ»(x)- p(x) =0, 


24 C33 233 + e33 dx 
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2 8033233 + 8e33 8 
€33 ( 129 Е o Я 
€33833 F €33 833233 + 8e33 8 
833 5 z 5 Ф› (х) + 
C33 (c33833 + езз )h h 


2 2 
t -2и Lon : uc ML. d h’ d (x) + 
2 2 (casgss F ei )л h dx 


6 
(свеззй? —c33e31h" ) ( ) 
€33 Сіз 
233 (свеззй? —c33e31h”) (casgas ce) d? 
+ is 7022 (х) = 
(c33833 ce)? C33 dx 
а? 
, giu (свеззй? — cex hi? (£u: z) 
x 
OT zx Q. 


54 2 
24 C33g33 + €33 
Задана разность потенциалов, механическая нагрузка равна нулю. Предполагается независимость попе- 
речного смещения от толщины и равенство нулю продольного перемещения и квадратичное распределение элек- 
трического потенциала по толщине (3): 
Из = (772 (х), ш = 0. (7) 


Выражение для поперечной силы: 
d d 
== —UZ h- —Ф р. 8 
О e 2 2) e (4 2 2) (8) 


C учетом равенства нулю нормальной компоненты вектора электрической индукции Hà лицевых IIOBepXHO- 
стях (Z = +h), уравнения для неизвестного прогиба UZ2(x) и распределения электрического потенциала (x) 


имеют вид: 
а? d’ 2 
=C44 —, UZ: (x) h-es — ®; (x) h—- p(x)-W phUZ; (x) = 0, (9) 
dx dx 
2 2 
02 (а) ви (1) -0. (10) 


Результаты исследования. Проведено сравнение результатов расчетов по предложенным прикладным 
теориям с расчетами колебаний пьезоэлемента (/ = 0,1 м, л = 0,01 м) на частоте равной 100 с! методом конечных 


элементов в АСЕГАМ [16]. 
В первой задаче, определяемой уравнениями (5) и (6), при задании равномерно распределенной нагрузки 


р(х) = 1000 Па:м и с граничными условиями: 
Zuz, (0)=0, О! k-0= 0, Ф (0) 2 ®2 (7) 2 0, UZ2 (1) = 0, Zuz, (1)=0 (11) 


получены следующие результаты, представленные на рис. 3. 4. 
Расчеты показали, что погрешность в определении вертикального смещения — 5,8 %, а для горизонтального 


смещения составляет 1,2 %. 


4,46 1Е-7 N 
Us, M 
2,788E-7 3x107 
2x107 
1,394E-7 3 
1x107 
0,000E+0 
0 


0,02 0,06 X M 
а) 6) 


Рис. 3. Вертикальное смещение в АСЕГАМ в первой задаче: 
а — распределение; 6 — график на верхней границе 
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Рис. 4. Горизонтальное смещение в АСЕГАМ в первой задаче: 
а — распределение; 6 — график на верхней границе 


Во второй задаче, определяемой уравнениями (8) и (9), при задании нулевой нагрузки p(x) = 0, разности по- 
тенциалов Vo = 100 В и граничными условиями: 


Qi k-0= 0, Ф (0) 20,05 (7) = Vo, UZ2 (I) 20 (12) 


получены следующие результаты, представленные Ha рис. 5, 6. 
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Рис. 5. Вертикальное смещение B ACELAN во второй задаче: 
а — распределение; 6 — график на верхней границе 
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Рис. 6. Электрический потенциал в АСЕГАМ во второй задаче: 
а — распределение; 6 — зависимость от продольной координаты в середине толщины 


Расчеты показали, что погрешность в определении вертикального смещения составляет 0,8 %, а для электри- 
ческого потенциала — менее | %. Следует отметить, что значения горизонтального смещения, рассчитанные в 
ACELAN, оказались на три порядка меньше максимального вертикального смещения, что говорит об адекватно- 
сти гипотезы (7). 

При задании механических нагрузок и разности потенциалов погрешность предложенного метода оказалась 
порядка 6 % для компонент перемещения и электрического потенциала. 
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Обсуждение и заключение. Как отмечалось в цитируемой литературе, одновременное использование изгиба 
и сдвига пьезоэлемента может значительно повысить его эффективность. Кроме этого, использование пористой 
керамики в силу различных зависимостей упругих модулей и пьезомодулей от процента пористости также 
улучшает выходные характеристики ПЭГ. 

В настоящей работе в силу линейной постановки в теории электроупругости удалось построить прикладную тео- 
рию расчета изгибно-сдвиговых колебаний пьезоэлемента в низкочастотной области, которая состоит из решения двух 
задач: в первой действуют механические нагрузки при нулевых потенциалах, а во второй — наоборот — нулевые 
механические нагрузки и задана разность потенциалов. На основе различных гипотез о распределении механического 
и электрического полей получены две краевые задачи для систем обыкновенных дифференциальных уравнений, ко- 
торые решались аналитически. Проведено сравнение результатов расчета смещений и электрического потенциала по 
предложенному методу и с помощью MKO, реализованного в пакете ACELAN. Эти расчеты подтвердили примени- 
мость предложенного метода, погрешность для которого в вычислении вышеуказанных характеристик составила 6 %. 
Такой точности достаточно для инженерных расчетов, поэтому предложенный метод может быть применен при IpO- 
ектировании пьезоэлектрических устройств, в том числе при сборе и накоплении энергии. Дальнейшее развитие этой 
прикладной теории будет направлено на охват более широкого частотного диапазона, включая первый изгибно-сдви- 
говый резонанс. 
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Аннотация 

Введение. Опубликованные исследования жесткости консолей под нагрузкой фокусируются на вопросах их де- 
формации и разрушения. Описаны расчеты момента инерции — принципиально важной характеристики проч- 
ности стержня. Однако не решена проблема значительных затрат времени для таких вычислений. Представлен- 
ное исследование восполняет данный пробел. Цель работы — описание нового быстрого метода аналитического 
расчета распределения напряжения сдвига в сечении консоли, соответствующего действию внешней приложен- 
ной силы. Впервые в таком контексте рассматриваются касательные напряжения и приводятся примеры расчета 
момента инерции для двух нестандартных сечений консоли. 

Материалы и методы. Для создания нового метода консоль представили как пачку пластинок, ориентирован- 
ных параллельно вектору внешней силы. Исходные расчеты строили по схеме консольной балки с выделенной 
пластинкой. Деформацию стержневых элементов моделировали с учетом действия однородного поля напряже- 
ния сдвига в сечении пластинки. Для обоснования упрощенного расчета момента инерции сечений задействовали 
схемы квадрата, эллипса, треугольника, шестиугольника, шестиконечной звезды и фигурного креста. Использо- 
вали аналитические и математические методы исследования, в частности теорему Гюйгенса-Штейнера. 
Результаты исследования. Создан быстрый универсальный метод вычислений момента инерции поперечного 
сечения консоли под нагрузкой. Его отличие — отказ от расчетов для каждого сечения с учетом формы и других 
особенностей. При любой форме сечения балка представляется как пачка бесконечно тонких пластинок, мо- 
менты их инерции интегрируются, и используется известное решение для прогиба тонкой пластинки. Метод поз- 
воляет однозначно показать распределение касательных напряжений на торце консоли, обеспечивающих задан- 
ный прогиб, причем впервые для таких решений используются касательные напряжения. Получены их профили 
в зависимости от направления внешней приложенной силы. Впервые выведены формулы для моментов инерции 
сложных сечений — шестиконечной звезды и фигурного креста. Каждое сечение соотнесено с кривой распреде- 
ления напряжения и его максимальным значением. Эти данные визуализированы в виде диаграмм. Установлено, 
что момент инерции и жесткость консоли не меняются при повороте внешней приложенной силы на 30° для 
сечения в виде звезды и на 45° — для квадрата и фигурного креста. В общем случае поле касательных зависит 
от геометрической формы и от ориентации сечения относительно внешней приложенной силы. 

Обсуждение и заключение. Предложенный упрощенный подход к расчету момента инерции поперечных сече- 
ний консолей дает возможность однозначно определить поле касательных напряжений на торце, обеспечиваю- 
щее при заданном прогибе соответствующее значение внешней приложенной силы. Инженеры и механики могут 
использовать результаты представленной работы при расчетах и моделировании деформации стержневых эле- 
ментов конструкций. 


Ключевые слова: деформация стержня, момент инерции плоской фигуры, момент инерции сложных сечений, 
упругий прогиб консоли, распределение касательных напряжений 
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Abstract 

Introduction. Published studies on the rigidity of consoles under load focus on the issues of their deformation and 
destruction. Calculations of the inertia moment, fundamentally important characteristic of the strength of the rod, are 
described. However, the problem of significant time consumption for such calculations has not been solved. The presented 
study meets the lack. The objective of the work is to describe a new rapid method for analytical calculation of the shear 
stress distribution in the section of the console corresponding to the action of an external applied force. For the first time, 
tangential stresses are considered, and examples of calculating the inertia moment for two non-standard sections of the 
console are given in this context. 

Materials and Methods. To develop a new method, the console was presented as a pack of plates oriented parallel to the 
vector of external force. The source calculations were based on the scheme of a console beam with a dedicated plate. The 
deformation of the rod elements was modeled taking into account the effect of a uniform shear stress field in the plate 
section. To validate the simplified calculation of the inertia moment of the sections, schemes of a square, ellipse, triangle, 
hexagon, six-pointed star, and a figured cross were used. Analytical and mathematical research methods were applied, 
specifically, the Huygens-Steiner theorem. 

Results. A rapid valid method for calculating the inertia moment of the cross section of the console under loading has 
been developed. Its difference is the rejection of calculations for each section, taking into account the shape and other 
features. For any shape of the section, the beam is represented as a bundle of infinitely thin plates, their inertia moments 
are integrated, and a well-known solution for deflection of a thin plate is used. The method allows us to unambiguously 
show the distribution of tangential stresses at the end of the console, providing a given deflection, and tangential stresses 
are used for such solutions for the first time. Their profiles are obtained depending on the direction of the external applied 
force. Formulas for the inertia moments of complex sections — a six-pointed star and a figured cross — are derived for 
the first time. Each section is correlated with the stress distribution curve and its maximum value. This data is visualized 
in the form of diagrams. It is found that the inertia moment and the rigidity of the console do not change when the external 
applied force is rotated by 30° for a star-shaped section and by 45° for a square and a figured cross. In general, the tangent 
field depends on the geometry and on the orientation of the section relative to the external applied force. 

Discussion and Conclusion. The proposed simplified approach to calculating the inertia moment of the cross sections of 
the consoles makes it possible to uniquely determine the field of tangential stresses at the end, which provides the 
appropriate value of the external applied force for the given deflection. Engineers and mechanics can use the results of 
the presented work in the calculations and modeling of deformation of rod structural elements. 


Keywords: rod deformation, inertia moment of a flat figure, inertia moment of complex sections, elastic deflection of the 
console, shear stress distribution 
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Введение. Многие строительные конструкции содержат элементы в виде стержней, которые в процессе 
изготовления или эксплуатации испытывают упругие деформации [1]. Прочность стержня или балки на изгиб 
определяет несущую способность конструкции [2]. Способность балки к упругой деформации характеризуется 
жесткостью, определяемой как отношение нагрузки P к упругому прогибу балки Ae : Mx = P / àe [3]. Как правило, 
в лабораторных условиях жесткость проверяют на консольной балке. Один ее конец заделывают в жесткое 
основание, а на другой действует внешняя сила, направленная перпендикулярно оси балки [4]. Моделирование и 
расчеты характеристик деформации и разрушения стержней, как правило, связаны с решением задачи о прогибе 
консольной балки, или консоли, под действием внешней приложенной силы [5]. При этом нет публикаций о 
простых и универсальных методах определения моментов инерции сложных сечений относительно действия 
внешней приложенной силы. Решение этой задачи описано в представленной статье. 

Цель данной работы — создание универсального, быстрого метода расчета моментов инерции сложных сече- 
ний консоли под действием внешней приложенной силы. Новый подход позволяет аналитически узнавать рас- 
пределение напряжения сдвига в сечении, соответствующее действию внешней приложенной силы. Отметим, 
что ранее в таких вычислениях не принимались во внимание касательные напряжения. Кроме того, впервые при- 
водятся примеры расчета момента инерции для сложных фигурных сечений. 

Материалы и методы. Во многих работах по сопротивлению материалов, например в [6], приводится уни- 
версальная формула для расчета упругого прогиба консоли, Ac. Согласно этой формуле, жесткость консоли равна: 

P/d, = ЗЕ, / 13, (1) 
где Е — модуль Юнга; L — длина консоли; /. — момент инерции поперечного сечения балки относительно OCH х, 
проходящей через центр тяжести сечения перпендикулярно приложенной силе Р. 

Из уравнения (1) следует, что принципиально важной характеристикой консоли является момент инерции 
сечения /,, величина которого зависит от геометрической формы поперечного сечения балки и направления 
ocu x [7]. Следует подчеркнуть, что в уравнении (1) момент инерции /. относится к оси X, которая перпендику- 
лярна направлению внешней приложенной силы P. В частности, момент инерции прямоугольного сечения axb 
относительно оси симметрии х равен [8]: 

I, - ab? /12. (2) 

Здесь а — толщина консоли; b — ее ширина. Сила Р направлена параллельно стороне b прямоугольника. 

Подставляя (2) в (1), получим известное уравнение для прогиба консоли прямоугольной формы [9]: 


3 
p 
Ea\ b 


Поперечные сечения консольных балок, или консолей, бывают разными. Ha рис. | приведен простой пример 
консоли квадратного сечения под действием внешней силы P, направленной вдоль диагонали квадрата. 


Рис. 1. Схема консольной балки с выделенной пластинкой 


Моментом инерции сечения относительно оси х называется сумма, или интеграл, произведений элементарных 
площадок ds = dxdy на квадраты расстояний у площадок до ocu x: 7, = [у?ахау. [10]. Подынтегральная функция, 


в сущности, является моментом инерции элементарной площадки dxdy относительно OCH X. 

Консоль можно представить в виде пачки предельно тонких пластинок толщиной dx и длиной Г, 
ориентированных параллельно вектору силы Р. Под действием Р все пластинки прогибаются на одну и ту же 
величину Ac. При заданной ориентации сечения консоли отдельная пластинка не испытывает влияния 
деформации остального объема. Тогда момент инерции сечения консоли в целом будет определяться 
интегральной суммой моментов инерции сечений всех пластинок в пачке. 

Проекция пластинки на плоскость поперечного сечения консоли представляет собой прямоугольную полоску 
толщиной dx и полудлиной h (рис. 1). Момент инерции сечения отдельной пластинки можно рассматривать как 
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момент инерции консоли прямоугольного сечения 2йхах. По определению, к каждой пластинке в пачке приме- 
нимо уравнение типа (2), где а=ах и b=2h. Согласно этому выражению, момент инерции полоски равен 
АГ. = 2h7(x)dx / 3. Таким образом, момент инерции сечения консоли можно определить интегрированием момен- 
тов инерции не элементарных площадок, а элементарных полосок: 


pus Е fiha’ ax, (3) 


где х изменяется в пределах от A до В. 

Условие ориентации плоскости пластин параллельно вектору внешней приложенной силы важно, поскольку позво- 
ляет однозначно связать упругий прогиб консоли Àe с распределением касательных напряжений в поперечном сечении 
консоли. Все пластинки в пачке прогибаются на одну и Ty же величину Ae Согласно (1), 
dl, = 2I? (x)dx / 3.Значит, для пластинки толщиной dx потребуется элементарная сила dP = 31.Е41,/ [2 = 21ЕЙ(х)аХ / L^. Ta- 
кая сила соответствует действию однородного поля напряжения сдвига в сечении пластинки ds(x) = 2h(x)dx: 

т=аР/ ds = A,ER (х) / D. (4) 

Из (4) видно, что значение напряжения t в системе координат Xy не зависит от координаты y. 

Уравнение (4) удобно использовать при моделировании деформации стержневых элементов конструкции. 

Интеграл (3) определяет момент инерции сечения /, относительно оси x, проходящей через центр тяжести сече- 
ния. В случае несимметричных и сложных сечений удобно сначала найти момент инерции сечения или части сече- 
ния относительно оси, которая не проходит через центр тяжести сечения. Затем нужно перейти к моменту инерции 
сечения J, относительно оси, которая проходит через центр тяжести сечения. Известно, что момент инерции сечения 
повторяет свойства момента инерции твердого тела и подчиняется теореме Гюйгенса-Штейнера [11]. Момент инер- 
ции сечения ly относительно произвольной оси x равен сумме момента инерции этого сечения J, относительно оси, 
проходящей через центр тяжести сечения параллельно оси х, и произведения площади сечения 5 на квадрат рассто- 
яния а между осями: 7, = [. + a?S. Поэтому в общем случае можно записать: 


_ 2 [8 3 Dea 2 
==], A(x) ах+а?5 = I, +a°S. (5) 


Если ось х проходит через центр тяжести сечения, то расстояние а = 0 и уравнение (5) переходит в (3). 

Результаты исследования. Для создания простого, быстрого и универсального метода вычислений отка- 
жемся от расчетов для каждого сечения с учетом его формы и других особенностей. Именно такой подход впер- 
вые реализован в рамках данной научной работы. Каким бы сложным ни было сечение, достаточно использовать 
известное решение для прогиба тонкой пластинки, представить балку в виде пачки бесконечно тонких пластинок 
и проинтегрировать их моменты инерции. Кроме того, метод позволяет однозначно показать распределение 
касательных напряжений на торце консоли, обеспечивающих заданный прогиб. Подчеркнем, что касательные 
напряжения в таком контексте рассматриваются впервые. 

Упрощенный расчет момента инерции простых сечений. Для обоснования предлагаемого метода рассмот- 
рим известные сечения простой геометрической формы. Далее вместо выражения «момент инерции поперечного 
сечения консоли» используем термин «момент инерции». Будем считать, что внешняя приложенная сила всегда 
направлена перпендикулярно оси х, относительно которой определяется момент инерции сечения. 

Квадрат. Уравнение (2) получено при условии, что вектор силы, приложенной к торцу консоли, 
перпендикулярен стороне а прямоугольника. При b = а получим момент инерции квадратного сечения Г = 44/12. 

Найдем предлагаемым методом момент инерции квадрата относительно оси х, которая параллельна не сто- 
роне, а диагонали квадрата (рис. 2 а). 


т, МРа 
у 
30 
20 
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a 10 
0 
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Рис. 2. Схема сечения в виде квадрата: а — сечение, 6— распределение напряжения т вдоль оси X 
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Полудлина полоски на рис. | равна л = x. Из уравнения (3) находим: 
2 ра/ , a^ 


= — хах = —. 
3 -a / 43 12 


Видно что поворот оси x на 45° не меняет момент ине ции квадратного сечения. Следовательно при повороте 
, , 


c 


внешней приложенной силы на 45? ne меняется и жесткость консоли (1). 

Подставляя x в вместо й (4), получим распределение напряжения сдвига в сечении консоли квадратного сече- 
ния, соответствующее прогибу Àe: 

(x) =x BID. 

На рис. 2 б представлено распределение т (х) вдоль диагонали квадрата при Ae = 2 MM Ha = 5 мм. Согласно (1), дей- 
ствие напряжения сдвига T(x) на торце консоли из стали (E = 200 MPa) длиной L= 50 мм соответствует действию 
внешней приложенной силы P = 0,25Е).а\Г?З = 500 N. Максимальное напряжение Tmax = 40 MPa наблюдается на Bep- 
тикальной диагонали квадрата. При отклонении от диагонали резко до нуля уменьшается напряжение т. В результате 
в системе т(х) образуется острый пик. 

Напряжение т зависит только от переменной x, поэтому по графику на рис. 2 6 можно определить значение T 
в любой точке квадрата. 

Очевидно, что при ориентации силы перпендикулярно сторонам квадрата x = Р/а? = 20 MPa (рис. 2 6, пунктирная ли- 
ния). Как видно, распределение напряжения T(x) в сечении, соответствующее внешней силе P, существенно зависит от ее 
направления. 

Эллипс. На рис. 3 а приведена схема эллипса с полуосями a и b. Начало координат — в центре эллипса. 


т, МРа 


60 
40 
20 
4 -2 —1 0 1 X, MM 
a) 6) 


Рис. 3. Схема сечения в виде эллипса c полуосями a и b: 


а — эллипс; 6 — распределения т вдоль полуоси а. 1 — a = b, 2 — b —2a, a = 2,5 MM 


Из канонического уравнения эллипса [12] следует, что полоска, выделенная на рис. 3, имеет полудлину 
h = Ба? — x)? / а. Подставляя эту величину в (3) и интегрируя от —а до +a, получим момент инерции эллиптиче- 
ского сечения: 


а 3 ра 3/2 Р 
I, == | Ware. [ a? -x | dx = v | (6) 
3 J-a За J-a 4 


При a = = гвыведем момент инерции кругового сечения: лғ“ / 4, где г — радиус окружности. 

Подставляя B (4) значение h, соответствующее эллипсу, получим следующее распределение напряжения 
сдвига в сечении консоли эллиптической формы: 

(x)= b (1 — х2 / a?)E | І. 

На рис. 3 6 представлено распределение T(x) вдоль полуоси а = 2,5 мм при Ae = 2 MM, E = 200 MPa и L = 50 мм. 
При сравнении ero с распределением t(x) на рис. 2 6 для квадрата очевидно существенное влияние геометриче- 
ской формы сечения консоли на распределение напряжений. В случае эллиптического сечения на кривой T(x) 
отсутствует острый пик. Максимальное напряжение Tmax наблюдается вдоль полуоси b. С учетом требования 
Ле = Const можно говорить о быстром росте напряжения Tmax при увеличении отношения b/a. При увеличении по- 
луоси b в два раза Tmax увеличивается в 4 раза (рис. 3 6). 


Механика 


163 


https://vestnik-donstu.ru 


164 


Дерюгин Е.Е. Упрощенный расчет момента инерции поперечного сечения консоли под нагрузкой 


Треугольник. Рассмотрим сечение в виде равнобедренного треугольника (рис. 4). Ось х направлена вдоль высоты 
треугольника (рис. 4 а). На расстоянии x от основания треугольника полудлина полоски равна й = b(a — x) / 2a. 


Рис. 4. Схема сечения в виде равнобедренного треугольника: 
а — ось х, сила направлена параллельно основанию; б — сила направлена перпендикулярно основанию 


Интегрирование выражения (4) похот0 до а определит следующую величину момента инерции: 


bb pe 3 bea 
а 7 
Ta SEPT (0) 


Рассмотрим случай, когда ось х проходит через центр тяжести и параллельна основанию треугольника (рис. 4 6). 


Полоска на расстоянии x от левого угла треугольника имеет полудлину h = х\/4а? / b? —1/2 (рис. 4 6). Центр 
тяжести полоски расположен на расстоянии / от основания треугольника (рис. 4 б). Следовательно, согласно тео- 


реме Гюйгенса-Штейнера для поперечного сечения, момент инерции полоски относительно основания треуголь- 
ника равен: 


L= 


3 
dics ТУРО ЧЕ ge 
3 3p 


Интегрирование полученного выражения по переменной x oT —b/2 до b/2 определит момент инерции треуголь- 
ника относительно основания: 


І, =a°b/12. (8) 


Центр тяжести треугольника расположен на расстоянии 4/3 от основания. Согласно теореме Гюйгенса-Штей- 
нера [13], момент инерции треугольника относительно собственного центра тяжести будет меньше (8) на величину: 


2 3 
r-(£ g-2b 
3 18 


Следовательно, момент инерции треугольника относительно собственного центра тяжести равен: 
ab ab ab (9) 
12 18 36` 

Полученный результат в точности соответствует табличному значению момента инерции сечения OTHOCH- 
тельно оси через центр тяжести, параллельной основанию треугольника. 


На рис. 5 приведены распределения напряжения сдвига в треугольном сечении при A, = 2 MM для случаев, 
когда сила прогиба направлена вдоль основания (а) и вдоль высоты треугольника (6). 


где 5 = ab/2 — площадь треугольника. 


c 
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Puc. 5. Распределения т в треугольном сечении: а — сила направлена вдоль основания; 
6 — сила направлена вдоль высоты треугольника 
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При ориентации треугольника как Ha рис. 4 а напряжение T(x) плавно уменьшается OT Tmax на основании до 0 
на вершине (рис. 5 а). При ориентации треугольника как на рис. 4 6 распределение t(x) (рис. 5 6) подобно распре- 
делению для квадратного сечения, когда сила ориентирована вдоль диагонали (рис. 2 6). 

Правильный шестиугольник. На рис. 6 представлены две ориентации шестиугольника относительно OCH X: 
параллельно (a) и перпендикулярно (6) его диагонали. 


dx 


6) 


Puc. 6. Правильный шестиугольник: 
а — диагональ перпендикулярна вектору силы; Ó — диагональ параллельна вектору силы 


Сечение шестиугольника на рис. 6 а состоит из: 

— прямоугольника шириной а и высотой av3; 

— двух равнобедренных треугольников высотой 4/2 и с основанием 4V3. 

Найдем моменты инерции указанных частей сечения по (2) для прямоугольника и по (7) для треугольников. 
Для прямоугольника в уравнении (2) b = av3, поэтому его момент инерции равен L4 = N3a*/4. Для треугольных 
частей в уравнении (7) основание равно b = av3, а высота a/2. Следовательно, момент инерции шестиугольника: 

І, = I4 +I = N3a* (1441/16) = 5Үза*/16. (10) 


Если приложенная сила направлена вдоль диагонали шестиугольника (рис. 6 6), TO полудлина полоски равна 
a/2 — x/N3. Совершив соответствующие замены в (4), получим: 


2 a43/2 543a* 
І, = 2 -x/ 4B dx =~. 11 
31.5, 0 Дал г" (10) 


Сравнивая (10) и (11), убеждаемся, что поворот шестиугольника на 30° не влияет на его момент инерции от- 
носительно оси х. 

Рассмотрим случай ориентации шестиугольника как на рис. 6 а. При = 2 мм и E = 200 MPa в интервале 
от —a/2 до +а/2 напряжение t(x) будет постоянным и равно Tmax = 15 MPa (рис. 7 а). При тех же условиях и ори- 
ентации шестиугольника как на рис. 6 б распределение т (х) подобно распределению для квадратного сече- 
ния (рис. 2 б). Однако здесь Tmax = 20 MPa И Tmin = 5 MPa (рис. 7 6). 
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Рис. 7. Распределение T(x) в сечении правильного шестиугольника: 


а — диагональ перпендикулярна вектору силы; 6— диагональ параллельна вектору силы 


Примеры упрощенного расчета момента инерции сложных сечений 

Шестиконечная звезда. Рассчитаем предлагаемым методом момент инерции нестандартного сечения в виде 
правильной 6-конечной звезды со стороной а. Для определения момента инерции относительно малой диагонали 
звезды воспользуемся схемой на рис. 8 а. Выделим на схеме три зоны: зону І шириной 4/2, зону П для остальной 
части полфигуры и примыкающую к ней вспомогательную зону Ш в виде треугольника. 
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45 
2 


a) 6) 
Рис. 8. К расчету момента инерции шестиконечной звезды относительно: 
а — малой диагонали; 6 — большой диагонали 


Полудлина полоски в зоне I на расстоянии x от диагонали равна hy = N3(a — x). Используя выражение (3), для 
двух 30H I на рис. 6 а получим момент инерции: 


а/2 И 
Га = 23] T T = 5 44, 
a/2 16 


Момент инерции зоны П равен моменту инерции прямоугольника высотой av3 и шириной а без момента 
инерции треугольника Ш. Полудлина полоски в прямоугольнике равна hy = a 3/2. Из уравнения (3) находим 
момент инерции прямоугольника: 


3 
[аш = at. 


4 
Полудлина полоски в треугольнике равна Am = х3. Согласно (4), момент инерции треугольника равен: 


1 E 
"1 


Следовательно, момент инерции зоны П равен: 


7T43 д 
Ш = Паш = Ги = uw E 
Удвоенная сумма моментов инерции зон I и П определит момент инерции 6-конечной звезды: 
11-43 
= таа а“. (12) 


Для определения момента инерции звезды относительно большой диагонали следует воспользоваться схемой 
на рис. 8 6, где выделены две зоны. Полудлина полоски в зоне I равна hy = a + x/N3, в зоне II hy = a/2 — x/3. Из (3) 


вычисляются моменты инерции зон I (Ja = a^65N3/96) и II (Ian = a*V3/96). Можно убедиться, что удвоенная сумма 
моментов инерции зон Ги П в точности соответствует уравнению (12). Следовательно, момент инерции 6-конеч- 


I xl 


ной звезды относительно большой и малой диагоналей одинаков. 

На рис. 9 а показано распределение напряжения сдвига в сечении 6-угольной звезды, соответствующее 
À — 2 MM и ориентации внешней приложенной силы согласно рис. 8 а. Вдоль вертикальной оси звезды напряже- 
ние принимает максимальное значение Tmax = 37,5 MPa. По мере удаления от оси вправо т быстро падает до 


уровня т = 9,375 MPa и остается постоянным в интервале 2 < x < 4 мм. Затем т падает до нуля. При удалении от 
оси влево т меняется аналогично. 
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Puc. 9. Распределения TB 6-угольном сечении: а — сила направлена вдоль большой диагонали; 
6 — сила на правлена вдоль малой диагонали 
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При повороте звезды Ha 30° распределение напряжения т существенно меняется, приобретая очертания ла- 
сточкиного хвоста (рис. 9 6). Наблюдаются два пика со значением Tmax = 28,12 MPa. Картина распределения сим- 
метрична. Однако при удалении от оси симметрии напряжение сначала увеличивается от 12,5 до 28,12 МРа, а 
затем падает до 3,125 MPa. В связи с этим на зависимости T(x) наблюдаются два пика. Далее напряжения быстро 
уменьшается до нуля. 

Фигурный крест. Рассмотрим решение для момента инерции нестандартного сечения в виде фигурного 
креста (рис. 10), каждая сторона которого является четвертой частью окружности радиуса К. Расстояние между 
вершинами фигуры равно RN2. Найдем сначала момент инерции креста относительно большой диагонали (рис. 10 а). 


a 
\ 


a) 


Рис. 10. Фигурный крест: a — сила направлена вдоль диагонали креста; 
6 — сила направленна под углом 45° к диагонали креста 


Половина длины полоски й на рис. 10 а равна 


h-R-4R? -(R-xy.. 


Учитывая симметрию и используя уравнение (3), получим следующее значение момента инерции сечения 
фигурного креста: 


R 3 
ы) Е Jx(2R x)] dx - [4-5x/ 4]R*. (13) 


Определим момент инерции креста при повороте ero на 45? относительно OCH x (рис. 10 6). 

Фигурный крест вписывается в окружность радиуса К. Из рис. 10 6 видно, что вокруг креста выделяются че- 
тыре фигуры в виде овала с острыми углами. Для определения момента инерции креста достаточно из момента 
инерции круга вычесть моменты инерции четырех овалов. 

Из (6) следует, что момент инерции круга равен: 


Lo = 1844. 


Определим моменты инерции овалов относительно центра тяжести креста. 
Момент инерции овала справа равен моменту инерции овала слева. Полувысота полоски у овала слева: 


h = JR? -(R-x)?. 


Учтем уравнение (3), a также симметричность данных частей и их расположения. Интегрированием от 0 до А — R/ V2 
получим следующее значение момента инерции данных двух частей: 


3 e R(1-1/42 ) 7 4 
pa [ i= a=x7R) | ax 8212 (14) 
3 Jo 314 


Центры тяжести овалов сверху и снизу окружности расположены на расстоянии а = R/V2 от центра тяжести 
целой фигуры. Площадь этих двух частей равна 5 = А2(л — 2). 
По (4) вычислим момент инерции данной пары: 
I, 2 a?8 = R*[n/2-]]. 


Величина / в (3), согласно рис. 10 6, (овал внизу круга) равна л = 42a? —x? —а. Следовательно, для пары 
овалов сверху и снизу круга можно записать момент инерции относительно собственных центров тяжести: 


R/N2 à 
pom -Šef | i-e ву EU dx = R*[3n/4- 7/3]. (15) 
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Дерюгин Е.Е. Упрощенный расчет момента инерции поперечного сечения консоли под нагрузкой 


Определим результирующий момент инерции фигурного креста. Для этого из момента инерции (13) круга 


вычтем моменты инерции четырех овалов: 
То - Ги = 12-213 = [4 51/4] В“. 


Сравнивая полученный результат с результатом (13), убеждаемся, что момент инерции данного крестообраз- 
ного сечения консоли не меняется при повороте оси x на 45°. 

На puc. 11 а показано распределение напряжения сдвига в сечении фигурного креста, соответствующее 
À = 2 мм и действию внешней приложенной силы вдоль оси креста (рис. 10 a). Ha вертикальной оси креста наблю- 
дается острый пик напряжения, где Tmax = 80,0 MPa. 


т, МРа т, MPa F 


0 
-3 -2 -1 0 1] 2x, MM 


Рис. 11. Распределения тв сечении фигурного креста: а — сила направлена вдоль диагонали креста; 
6 — сила направлена под углом 45° к диагонали креста 


При повороте креста относительно приложенной силы на 45° распределение напряжения т существенно ме- 
няется. По мере отклонения оси симметрии напряжение сначала плавно увеличивается с т= 13,73 МРа до 
Tmax = 20,32 MPa. Затем напряжение т резко падает до нуля. Поэтому в распределении т наблюдается два симмет- 
ричных пика (рис. 11 6) на расстоянии 2 мм от оси симметрии. 

Обсуждение и заключение. Предложенный упрощенный метод расчета позволяет быстро узнать моменты 
инерции сложных поперечных сечений консоли. При этом однозначно определяется поле напряжения сдвига в 
сечении образца, соответствующее действию внешней приложенной силы. Кроме того, показано, что 
распределение напряжений в сечении качественно и количественно зависит от ориентации сечения относительно 
направления внешней приложенной силы. 

Для обоснования метода проведены вычисления моментов инерции не только для известных сечений простой 
геометрической формы (которые показали абсолютную идентичность расчетных и опубликованных в литературе 
результатов), но и для двух новых сложных сечений в виде правильной шестиконечной звезды и фигурного 
креста. Показано, что жесткость консоли не меняется при повороте внешней силы, приложенной 
перпендикулярно оси симметрии, на 30° для сечения в виде 6-конечной звезды, и на 45° — для квадрата и 
фигурного креста. Метод и полученные решения могут быть использованы инженерами и механиками при 
моделировании и расчетах прочности и жесткости стержневых элементов конструкций. 
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Аннотация 


Введение. Освоение полярных районов Мирового океана, необходимость решения различных задач, связанных 
с наличием большого числа замерзающих внутренних водоемов, ставят перед наукой новые проблемы. К их 
числу относится проблема изучения поведения ледяного покрова под воздействием на него различного вида 
нагрузок. Большой интерес представляет рассмотрение задач о действии на ледяной покров подвижной нагрузки. 
Подвижная нагрузка моделирует действие на лед движущихся транспортных средств. Однако в работах, посвя- 
щенных вышеуказанным задачам, рассматриваются случаи движения нагрузки по прямолинейной траектории. 
Целью данной работы является разработка метода исследования поведения ледяного покрова под действием 
нагрузки, перемещающейся произвольным образом. 

Материалы и методы. В статье предложен метод решения задачи о действии на ледяной покров водоема ко- 
нечной глубины движущейся по произвольной траектории силы. Задача сводится к решению системы двух диф- 
ференциальных уравнений. Первое из них моделирует поведение ледяного покрова и является уравнением коле- 
баний вязкоупругой пластины. Второе — моделирует поведение жидкости, находящейся в состоянии потенци- 
ального течения, и является уравнением Лапласа. Для решения системы дифференциальных уравнений приме- 
нялись интегральные преобразования по временной и пространственным переменным. Полученное в результате 
решение выражалось через повторный интеграл, для вычисления которого применялись численные методы. 
Результаты исследования. В результате реализации предложенного метода получено решение задачи о движе- 
нии сосредоточенной силы по ледяному покрову по произвольному закону. При этом произведены исследования 
характера поведения перемещений и напряжений в ледяном покрове в зависимости от скорости и ускорения дви- 
жения вертикальной нагрузки, глубины водоема и вязкоупругих свойств льда. Кроме того, рассчитано распреде- 
ление вектора скорости частиц жидкости по глубине водоема. 

Обсуждение и заключение. Предложенный метод является весьма эффективным для решения задач о подвиж- 
ных нагрузках, действующих на ледяной покров водоема конечной глубины. Он позволяет решать задачи о дей- 
ствии нагрузки, движущейся по ледяному покрову по сложной траектории. Полученные результаты могут быть 
использованы для расчета напряжения и перемещений ледового покрова при прокладке ледовых дорог или стро- 
ительстве аэродромов на льду. 


Ключевые слова: бесконечный ледяной покров, движущаяся нагрузка, произвольная траектория, переменная скорость 
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Abstract 

Introduction. The development of the polar areas of the World Ocean and the need to solve various problems associated 
with a large number of freezing inland water bodies issue new challenges for science. These challenges include the 
problem of studying the behavior of ice cover when exposed to various types of loads. Of great interest is the consideration 
of problems about the action of a moving load on the ice cover. A moving load simulates the effect of moving vehicles 
on ice. However, in papers devoted to the above problems, cases of load movement along a straight-line trajectory are 
considered. The objective of this research is to develop a method for studying the behavior of ice cover under the action 
of a load moving arbitrarily. 

Materials and Methods. The article proposes a method for solving the problem of the action of a force moving along an 
arbitrary trajectory on the ice cover of a reservoir of finite depth. The problem amounts to solving a system of two 
differential equations. The first of them models the behavior of the ice cover, and it is the equation of vibrations of a 
viscoelastic plate. The second equation simulates the behavior of fluid in a state of potential flow, and it is Laplace's 
equation. To solve the system of differential equations, integral transformations in time, space and variables were used. 
The resulting solution was expressed through an iterated integral, which was calculated using numerical methods. 
Results. The development and implementation of the method resulted in solving the problem of the movement of a 
concentrated force along an ice cover according to an arbitrary law. At the same time, studies were carried out on the 
behavior of displacements and stresses in the ice cover depending on the speed and acceleration of the movement of the 
vertical load, on the depth of the reservoir, and on the viscoelastic properties of ice. In addition, the distribution of the 
velocity vector of fluid particles along the depth of the reservoir was calculated. 

Discussion and Conclusion. The proposed method is very effective for solving problems of moving loads acting on the 
ice cover of a reservoir of finite depth. It provides solving problems about the action of a load moving along an ice cover 
along a complex trajectory. The results obtained can be used to calculate the stress and displacement of the ice cover 
during the laying of ice roads or the construction of airfields on the ice. 


Keywords: infinite ice cover, moving load, arbitrary trajectory, variable speed 
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Введение. Освоение полярных районов Мирового океана и наличие большого количества замерзающих BHYT- 
ренних водоемов приводят к необходимости изучения полей перемещений и напряжений ледяного покрова, обу- 
словленных действием различного вида нагрузок. Решению этих проблем посвящено большое количество работ 
отечественных и зарубежных ученых. Ранее было установлено, что механические свойства льда зависят от соле- 
ности воды и его температуры. Большое внимание уделялось разработке численных моделей льда, которые до- 
статочно точно отражали взаимодействие льда и идеальной несжимаемой жидкости. В работах [1, 2] для этого 
применялся метод гидродинамики сглаженных частиц, в [3, 4] — метод дискретных элементов. В статье [5] лед 
моделируется упругой пластиной, лежащей на поверхности стратифицированной жидкости. Модели, допускаю- 
щие наличие трещин, рассматривались в работах [6, 7]. Модели льда, усиленного армирующими элементами, 
представлены в работах [8, 9]. 

При этом в некоторых работах ледяной покров рассматривается как упругая пластина, лежащая на поверхно- 
сти водоема [10, 11]. В то же время в [12] на основе приведенных исследований делается вывод о том, что в 
некоторых случаях свойства льда лучше всего описывает реологическая модель Кельвина-Фойгта с одним пара- 
метром (временем затухания). Поэтому многие исследователи при моделировании ледового покрова применяют 
вязкоупругую пластину [13]. В [14] для описания свойств льда применялись нелинейные модели. 
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В некоторых работах рассматривалось действие Ha ледяной покров подвижной нагрузки. В [15] исследовалось 
действие подвижной нагрузки на ледяной покров в замороженном канале, в [16] рассматривалось действие на 
ледяной покров нагрузки с импульсивным характером движения. Работа [17] посвящена изучению нагрузки, дви- 
жущейся по замерзшему руслу. При этом исследовалось прямолинейное движение нагрузки [18]. Однако в ре- 
альных условиях часто приходится иметь дело с нагрузкой, перемещающейся более сложным образом. Поэтому 
целью данной работы явилась разработка метода решения задач о действии нагрузки, движущейся по ледяному 
покрову по сложной траектории. Это позволит более точно исследовать действие на лед транспортных средств, 
движущихся сложным образом. 

Данная работа является продолжением исследований, связанных с рассмотрением задач о действии подвиж- 
ной нагрузки на различные объекты, результаты которых представлены в работах [19, 20]. 

Материалы и методы. Постановка задачи. Рассматривается водоем конечной глубины с бесконечным ледя- 
ным покровом (бесконечная пластина), который подвержен действию вертикальной силы, движущейся произ- 
вольным образом — импульсивно. Предполагается, что жидкость водоема несжимаема и совершает потенциаль- 
ное движение. 

Задача сводится к системе дифференциальных уравнений [15]: 


а = 22275), 
AF =0, 
где W(x, y, — прогиб ледяного покрова; E и и — соответственно модуль Юнга и коэффициент Пуассона льда; 
D-ERB/I2( — p?) — цилиндрическая жесткость на изгиб; h — толщина ледяного покрова; to — время релаксации де- 
формаций; Aj = (022 + 67); A= 22 + 6,7 + 02; рли рь — соответственно плотность льда и воды; с? = ps VD; К = ps 2/0; 
b=p,/D; Q(x, y ,t) — действующая на поверхности льда нагрузка; F(x, у, Z,  — потенциал скорости. 
При краевых условиях при 2=0 (граница лед-вода): 
ôW = 0,Е. 


На дне водоема при z = —H: 

0,F =0. 

Кроме того, предполагалось, что ледяной покров и жидкость в водоеме в начальный момент времени находи- 
лись в состоянии покоя. Нагрузка представляла собой сосредоточенную единичную силу (величиной в один 
Ньютон) Q(x, y, t), которая перемещалась произвольно по незамкнутой кривой y произвольной формы. Полага- 
лось, что перемещение силы задавалось в виде О = Q(s(t)), где s — дуговая координата, отсчитываемая OT HEKO- 

Е T x 2 хо (t) 
торой фиксированной точки траектории y. Траектория движения задавалась параметрически в виде ; 
y 2 У, (1) 
где t — время. 
Подвижная нагрузка аппроксимировалась выражением 


O(x,y,t) ==? exp( e’ ((« Xo (0) t (y — уо (y )/. 


где = — числовой параметр. 
После применения интегрального преобразования Фурье по переменным х и у, интегрального преобразова- 
ния Лапласа по f, были получены формулы для вычисления неизвестных функций Wu F: 


1 t oo | А 1 
УР = 22-р /4& J R -(a-d) -(a+d)t 
W (x».)7 275 JJ» е о (р (r-3)- [e —e Japas, 


нон 


160 _ e C | рат, 


t oo 
1 2 2 
(х) | | pte in а): 
0 


0 
С + d)e (ray sue а) ] арат, 


R? (т) -8? +р?, 8 = xo (т) —х, В = yo(x)- у, 
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5 


i 1/2 Top 
Y - [sop –4(с2р? + bpcth( pH))(p* +k) | a= 2(c?p рН) 


а= l | 
2(c*p +beth( pH)) 

Используя известные соотношения из теории тонких пластин и теории потенциального течения идеальной 
жидкости, можно получить соотношения для вычисления перемещений и напряжений в ледовом покрове, а также 
компонент вектора скорости частиц жидкости. 

При вычислении несобственного интеграла с помощью численных методов использовалось приближенное 


oo A 
COOTHOIIICHHC | f(p)dp -f f(p)dp, в котором величина А выбиралась настолько большой, чтобы оценка 
0 0 


oo 
ошибки | f(p)dp не превышала установленного значения. 
А 


Так, для величины прогиба льда 


t œ tA t œ 
1 1 1 
W (x,y,t)= т |^ (т, р)арат = z5]J ^ (spi | [us (т, р) рат, 


данная оценка имеет вид: 


t oo 
1 2 АЕ? 214,2 
—— | | Uo (т, р) арат|  ——— С, 
sab | Ye p) Ki xDy(A)- 


y2 
v(4) = [84 -4(e? 4?  bAeth( AH )(A* +k) |. 
Подобные оценки можно получить и для остальных вычисляемых величин. Эти оценки использовались при 
определении величины A, 
2 
2 АЕ A2 [4=? 


В проведенных расчетах величина А выбиралась такой, чтобы оценка арта) 
TY 


не превышала 


1 t о 
—— Оз (т, p)dpdt| 0,001. 
551,1, of p) á 


При вычислении повторного интеграла использовались квадратурная формула Симпсона (по переменной T) и 
квадратурная формула Чебышева с равными весами для двух узлов (по переменной р). Аналогичным образом 
вычислялись и остальные величины. 

Результаты исследования. Разработан метод решения задач о действии нагрузки, движущейся по ледяному 
покрову водоема, наполненного идеальной жидкостью, по сложной траектории с переменной скоростью. Приме- 
няя данный метод, проведены расчеты, которые показали степень влияния различных параметров на деформацию 
ледяного покрова. 

Изложенный метод не накладывает ограничений на форму траектории движения сосредоточенной силы. При 
расчетах рассматривался частный случай траектории, состоящей из дуг окружностей (рис. 1). Красной точкой 
указано положение сосредоточенной силы в рассматриваемый момент времени и направление движения силы. 


У, Mo 


10 | 


66 


0 | 4 | 8 i 12 | 16 Xm 


Рис. 1. Траектория движения сосредоточенной силы 
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При расчетах принимались следующие значения используемых параметров: толщина ледяного покрова 
h = 0,25 м, модуль Юнга E = 500 000 000 Н/м?, коэффициент Пуассона льда u = 1/3, плотность льда р = 900 кг/м?, 
плотность воды р = 1 000 кг/м?, e = 5. Ниже представлены результаты проведенных расчетов. 

На рис. 2 приведено изменение прогиба ледяного покрова при скорости движения силы у = 2,5 м/с, касатель- 
ном ускорении у, = 1 m/c’, глубине водоема H = 25 м и времени релаксации To= 1 с. 

Закон движения силы по траектории принимался в виде: 

з= 4101 ast ay. 
Коэффициенты ai, a2, аз подбирались таким образом, чтобы сила, находясь B одной и той же точке траектории, 


имела необходимые скорость и касательное ускорение. 


Рис. 2. Изменение прогиба ледяного покрова 


При других значениях указанных параметров качественный характер распределения прогиба ледяного по- 


крова оставался практически неизменным. 


ZM 


Y,M 
Рис. 3. Движение жидкости, обусловленное действием подвижной нагрузки на ледяной покров 


Движение жидкости, обусловленное действием подвижной нагрузки при тех же значениях скорости, каса- 
тельного ускорения движения нагрузки, времени релаксации и глубины водоема, представлено на рис. 3 (изоб- 
ражено распределение вектора скоростей частиц жидкости). 


W, м:107 W м:107 


-9,0 -7,0 
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Рис. 4. Изменение величины максимального прогиба ледяного покрова в зависимости: 
а — от величины касательного ускорения; 6 — от глубины водоема 
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Влияние на максимальный прогиб ледяного покрова касательного ускорения движения силы представлено на 
рис. 4 а. В этом случае скорость движения силы была равной у = 17,5 м/с, а время релаксации то = | с. 

На рис. 4 6 изображен график зависимости величины максимального прогиба ледяного покрова W от глубины 
водоема Н. Скорость движения нагрузки при этом равнялась у = 17,5 м/с, касательное ускорение у; = 1 м/с? и 
время релаксации To= 1 с. 


FL gm 


-1,0 } 


1.4. 


J 
0 10 20 30 Г, м/с 
Рис. 5. Изменение максимального прогиба ледяного покрова в зависимости от скорости движения силы 


Зависимость максимального прогиба ледяного покрова от скорости движения силы представлена графически 
на рис. 5. Глубина водоема при этом принималась равной 25 M, а величина касательного ускорения — wr = | м/с2. 
Сплошной линией изображена зависимость, соответствующая времени релаксации т, = 1 с, пунктирная линия 
соответствует времени релаксации т, = 10 с. 

Обсуждение и заключение. Исследовано влияние глубины водоема на максимальный прогиб льда. Получена 
картина прогиба ледового покрова, обусловленная действием сосредоточенной силы, движущейся по сложной 
траектории с переменной скоростью. Расчеты показали, что с увеличением глубины водоема максимальный про- 
гиб ледяного покрова уменьшается (рис. 2). При этом заметная зависимость прогиба ледяного покрова от глу- 
бины водоема Н имеет место лишь для Н < 25 м. При больших глубинах величина максимальных прогибов ста- 
билизируется около некоторого постоянного значения и практически не изменяется. Таким образом, если 
Н> 25 м, то при расчетах глубину водоемов можно считать бесконечной. 

Увеличение касательного ускорения приводит к увеличению прогиба ледяного покрова. Причем зависимость 
прогиба от касательного ускорения очень близка к линейной зависимости (рис. 4). 

При малых временах релаксации To скорость движения нагрузки заметно влияет на величину прогиба льда. При 
больших временах влияние скорости движения нагрузки на прогиб ледяного покрова заметно уменьшается (рис. 5). 

Для изучения состояния жидкости водоема определено распределение вектора скорости движения частиц 
жидкости, обусловленное действием подвижной силы по льду (рис. 3). 

Разработанный метод решения задач и результаты, полученные с его помощью, могут быть использованы при 
строительстве ледовых дорог, проектировании и строительстве взлетно-посадочных полос на льду. 
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Аннотация 


Введение. Компьютерное моделирование позволяет инженерам принимать обоснованные проектные решения за 
счет точной оценки тепловых характеристик объектов проектирования. Актуальным направлением научных ис- 
следований и разработок является реализация технологии цифровых двойников в процессе проектирования тех- 
нических объектов. Для этого необходимо разрабатывать компьютерные модели, точность которых соответ- 
ствует требованиям, предъявляемым к цифровым двойникам. Однако в научной литературе недостаточно ши- 
роко представлены результаты исследований, направленных на реализацию технологии цифровых двойников в 
процессе проектирования. В основном рассматриваются общие вопросы, связанные с применением цифровых 
двойников в различных отраслях промышленности. Поэтому целью данного исследования явилась разработка 
цифровой модели и сравнительный анализ точности расчетов тепловых характеристик объекта проектирования. 
Материалы и методы. В качестве основного инструмента для проведения исследования выступает предложен- 
ная авторами методика разработки компьютерной модели тепловых характеристик для реализации технологии 
цифровых двойников. Численное решение реализовано путем построения тепловой модели для расчета темпера- 
турного поля на основе метода конечных элементов в системе инженерного анализа «Ansys» от компании «Ansys 
Inc» (США). Для аналитического решения применяется разработанная на основе метода пространства состояний 
компьютерная модель тепловых характеристик, реализованная в модуле «Ansys Twin Builder». Модель простран- 
ства состояний приводится в соответствие с поведением исходной тепловой модели путем приближения переда- 
точной функции к пошаговому отклику тепловой нагрузки с применением метода векторной аппроксимации во 
временной области. Верификация построенной аналитической модели выполнялась в системе инженерных рас- 
четов «MATLAB» от компании «The MathWorks» (США). Исследования проводились для станка модели 400V 
производства предприятия ООО «НПО «Станкостроение» г. Стерлитамак (Россия). 

Результаты исследования. Разработана цифровая модель, позволяющая с высокой точностью выполнить расчет теп- 
ловых характеристик объекта проектирования. Результаты сравнительного анализа показывают высокую степень со- 
ответствия значений тепловых характеристик, полученных с помощью предложенной цифровой модели, результатам 
численного моделирования. Максимальная погрешность расчета тепловых характеристик не превышает 0,1 °С. 
Обсуждение и заключение. Компьютерное моделирование, сочетающее численные методы расчета и научный 
подход, основанный на технологии цифровых двойников, позволяют получить результат максимально прибли- 
женный к результатам экспериментов. Предложенная в исследовании цифровая модель является эффективным 
решением, поскольку позволяет выполнить расчеты для оценки тепловых характеристик в режиме реального 
времени, что является одним из важнейших требований при реализации технологии цифровых двойников. 


Ключевые слова: цифровой двойник, сложный технический объект, компьютерное моделирование, 
автоматизированное проектирование, температурное поле, тепловые характеристики 
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Abstract 

Introduction. Computer modeling allows engineers to make valid design decisions by accurately assessing the thermal 
characteristics of design objects. The implementation of digital twin technology in the process of designing technical 
facilities is the current direction of scientific research and development. To do this, it is necessary to develop computer 
models whose accuracy meets the requirements for digital twins. However, the scientific literature does not widely present 
the results of research aimed at implementing digital twin technology in the design process. The general issues related to 
the use of digital twins in various industries are mainly considered. Therefore, the objective of this study was the 
development of a digital model and a comparative analysis of the accuracy of calculations of thermal characteristics of 
the design object. 

Materials and Methods. The main tool for conducting the research was the methodology proposed by the authors for 
developing a computer model of thermal characteristics for the implementation of digital twin technology. The numerical 
solution was implemented through constructing a thermal model for calculating the temperature field based on the finite 
element method in the ANSYS engineering analysis system from ANSYS, Inc. (USA). For the analytical solution, a 
computer model of thermal characteristics developed on the basis of the state-space method, implemented in the ANSYS 
Twin Builder module, was used. The state-space model was matched to the behavior of the original thermal model through 
approximating the transfer function to the stepwise response of the thermal load using the time domain vector 
approximation method. Verification of the constructed analytical model was carried out in the engineering calculation 
system MATLAB from the MathWorks company (USA). The research was carried out for a 400V machine model 
manufactured by NPO “Stankostroenie” LLC, Sterlitamak (Russia). 

Results. The developed digital model makes it possible to calculate the thermal characteristics of the design object with 
high accuracy. The results of the comparative analysis showed a high degree of correspondence between the values of 
thermal characteristics obtained using the proposed digital model and the results of numerical simulation. The maximum 
error in calculating thermal characteristics did not exceed 0.1°C. 

Discussion and Conclusion. Computer modeling that combines numerical calculation methods and a scientific approach based 
on digital twin technology, provides obtaining the result as close as possible to the results of experiments. The digital model 
proposed in the study is an effective solution, since it provides performing calculations to evaluate thermal characteristics in real 
time, which is one of the most important requirements for the implementation of digital twin technology. 


Keywords: digital twin, complex technical object, computer modeling, computer-aided design, temperature field, 
thermal characteristics 
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Введение. Компьютерное моделирование традиционно является эффективным инструментом для решения 
тепловых проблем на ранней стадии проектирования сложных технических объектов. Однако решение про- 
блемы требует точной оценки тепловых характеристик объекта проектирования для снижения негативных эф- 
фектов, вызванных повышением температуры [1]. При этом одним из эффективных инструментов для предва- 
рительной оценки тепловых характеристик является имитационное моделирование в системе инженерного 
анализа на основе передовых цифровых решений и разработок. Например, в работе [2] авторов Jianying Xiao и 
Kaiguo Fan представлен разработанный цифровой двойник для определения тепловых характеристик. Принцип 
работы двойника заключается в моделировании теплового поведения объекта путем отображения и коррекции 
тепловых граничных условий. Результаты эксперимента показали высокую точность модели (более 95 96), что 
существенно для повышения точности моделирования тепловых характеристик и тепловой оптимизации. По- 
этому актуальным направлением научных исследований и разработок в области моделирования является при- 
менение систем искусственного интеллекта [3] и цифровых двойников [4]. Например, в работе [3] авторов 
Наогап У1 и др. для повышения точности анализа тепловых характеристик предложена интерактивная модель 
коррекции тепловых граничных условий на основе нейронной сети. Результаты экспериментов показывают, 
что точность расчета температурного поля превышает 98 %, а точность прогнозирования тепловой деформа- 
ции — 96 %, что эффективно повышает точность моделирования. Кроме того, в работе [5] авторов Kypra- 
нова Н.В. и др. отмечено, что цифровые двойники часто применяются с целью совершенствования физических 
прототипов сложных технических объектов, поскольку позволяют не только обеспечить информационное со- 
провождение процесса проектирования, но и принимать эффективные проектные решения на основе наработок 
в области искусственного интеллекта. 

Одной из характерных особенностей технологии цифровых двойников является то, что для имитационного мо- 
делирования часто применяются компьютерные модели пониженного порядка [6]. Поэтому разработка компьютер- 
ных моделей является одним из базовых условий в реализации технологии цифровых двойников [7]. Понижение 
порядка моделей — эффективный и понятный с математической точки зрения подход для преодоления ограничений 
времени выполнения многомерных имитационных моделей. Например, в работе [8] авторов Мирзаев Д.А. и др. 
предложена простая аналитическая модель тепловых полей для разработки цифровых двойников процесса про- 
мышленной дуговой сварки. Авторы Бордачев Е.В. и Лапшин В.П. в своей работе [9] представили результаты ма- 
тематического моделирования температуры в зоне контакта инструмента и изделия при токарной обработке метал- 
лов. Такой подход позволяет получить точную оценку тепловых характеристик, соответствующую результатами 
численных экспериментов в режиме реального времени. В работе [10] авторы Schröder С. и Matthias У. представили 
модель пониженного порядка и предложили новую процедуру балансировки модели, основанную на преобразова- 
нии сдвига состояния. В заключении представили результаты сравнительного анализа с результатами из источников 
литературы, полученных с помощью серии численных экспериментов. 

Применение линейной и инвариантной по времени компьютерной модели пониженного порядка позво- 
ляет обеспечить быстрое моделирование при сохранении высокой точности расчетов [11]. При разработке 
компьютерной модели выполняется аппроксимация [12] передаточной функции для приближения модели 
пространства состояний к пошаговому отклику исходной тепловой модели [13]. Поскольку пошаговый от- 
клик тепловой нагрузки получен из исходной тепловой модели, то цифровая модель должна предоставить 
те же значения тепловых характеристик. 

Однако несмотря на то, что в последнее время наблюдается рост интереса к цифровому двойнику, в научной 
литературе недостаточно широко представлены результаты исследований, связанных с реализацией цифровых 
решений в процессе проектирования технических объектов. На основе систематического обзора литературы и 
тематического анализа публикаций по цифровым двойникам выявлен один из ключевых пробелов в знаниях, 
связанный с разработкой математического, программного и методического обеспечения высокоточных компью- 
терных моделей в рамках реализации технологии цифровых двойников. 

В этой связи целью исследования являлась разработка компьютерной модели пониженного порядка и анализ 
эффективности её применения в составе цифровой модели для точного расчета тепловых характеристик сложных 
технических объектов проектирования. 

Для достижения поставленной цели требовалось построить тепловую модель и выполнить расчет темпера- 
турного поля объекта проектирования, сгенерировать независимые пошаговые отклики тепловой нагрузки с по- 
мощью разработанного программного сценария [14], реализовать цифровую модель для расчета тепловых харак- 
теристик, определить погрешность расчетов, полученных с помощью численного и аналитического решений и 
провести сравнительный анализ полученных результатов моделирования. 


Позевалкин В.В. и др. Реализация цифровой модели тепловых характеристик на основе температурного поля 


Материалы и методы. Построение температурного поля объекта моделирования выполнялось для однород- 
ного изотропного тела на основе уравнения нестационарной теплопроводности: 


OT[0t = a? AT +9, с, (1) 


где T — температура (°С); t — время (c) a -4/A/(p-c) — коэффициент температуропроводности (м/с); 
X— коэффициент теплопроводности (Вт/м-°С); с — удельная теплоемкость (Дж/кг:°С); p— плотность материала (кг/м?); 
А — оператор Лапласа; д, — объемная мощность тепловыделения (Br/M?). 

Тепловой поток в процессе теплопередачи принимался равным количеству теплоты, переносимой через про- 
извольную поверхность площадью 5 в единицу времени f, что выражается следующим уравнением: 


js. | f аза, (2) 


где О — тепловой поток (Вт); qn — плотность теплового потока (Br/w?); S — площадь поверхности (м2). 

Плотность теплового потока при теплоотдаче определялась по формуле: 

qe =a-S-(Ts -Т.), (3) 

где a — коэффициент теплоотдачи (Вт/(м?-°С)); Ts — температура поверхности (°C); To — температура окру- 
жающей среды (°С). 

В этой связи, для расчета температурного поля объекта моделирования по формуле (1), задавались тепловые 
потоки (2) и (3), которые определяли количество теплоты, проходящее через поверхность в единицу времени. 

В качестве объекта моделирования был выбран компонент (рис. 1) металлорежущего станка модели «400V» 
производства предприятия ООО «НПО «Станкостроение» (Россия, г. Стерлитамак) в виде сверлильной головки. 


Рис. 1. Геометрическая модель объекта в «Ansys Design Modeler»: 
1 — корпус; 2 — электродвигатель; 3 — шпиндельный узел; О — тепловые потоки; 
qv — объемная мощность тепловыделений; qn — плотности тепловых потоков; 
a — коэффициент теплоотдачи; Ti-T4 — температурные датчики 


Поскольку для расчета температурного поля учитывалось количество теплоты, выделяемое в основном элек- 
тродвигателем (электромагнитные потери) и шпиндельным узлом (механические потери на трение), то в качестве 
внутренних источников тепла принимались электродвигатель qyi и шпиндельный узел 4,2, вблизи которых зада- 
вались соответствующие тепловые потоки О! и Q». Плотности тепловых потоков qni, 92 назначались для поверх- 
ностей, расположенных вблизи электродвигателя и шпиндельного узла соответственно. 

Конвекция задавалась коэффициентом теплоотдачи а с учетом условий теплообмена (свободная конвекция B 
воздухе). Поскольку станок контактирует с газообразной средой (воздух), то количество теплоты, отданное 
нагретой поверхностью окружающей среде в единицу времени f, прямо пропорционально разности температуры 
поверхности Ts и среды To в зависимости от 5 площади теплоотдающей поверхности (3). 

При построении тепловой модели объекта (твердого тела), состоящего из однородного материала (конструк- 
ционная сталь) с постоянными теплофизическими свойствами и наличием внутренних источников тепла, назна- 
чались следующие начальные и граничные условия: 

— начальные условия учитывали фиксацию постоянной температуры по всей поверхности объекта моделиро- 
вания (t= 0: T= То = const); 

— граничные условия второго рода задавались тепловыми потоками электродвигателя (Q1), шпиндельного узла (Q2), 
плотностью теплового потока (q,1) от электродвигателя к передней стенке корпуса и (4,2) к внутренним поверхностям; 

— граничные условия третьего рода задавались коэффициентом теплоотдачи (Q) для поверхностей, располо- 
женных внутри корпуса сверлильной головки; 

— граничные условия четвертого рода для контактных соединений поверхностей учитывали идеальный теп- 
ловой контакт и отсутствие термического сопротивления: 


Ta =i, +4,. (4) 
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Тепловые потоки (Qi, О»), плотности тепловых потоков (qni, 42), как и объемная мощность тепловыделения 
(qvi, 9), назначались в соответствии с известными рекомендациями по металлорежущим станкам [1]. Граничные 
условия, как и тепловые потоки, задавались для наружных поверхностей. Поэтому учитывалось, что взаимосвязи 
тепловых полей присутствовали только между внешними поверхностями. 

Дифференциальное уравнение (1), совместно с начальными и граничными условиями второго (2), третьего (3) 
и четвертого (4) рода, представляет собой математическую формулировку поставленной задачи. Поставленная 
задача решалась с помощью методов численного и аналитического моделирования. 

Численное решение выполнялось на основе метода конечных элементов в системе инженерного анализа 
«Ansys», которая разрабатывается американской компанией «Ansys Inc» (США) и поставляется фирмой «МЦД», 
авторизованным дистрибьютором «Ansys» в России. Геометрическая модель (рис. 1) объекта импортировалась в 
проект «Ansys Workbench» с добавлением блока анализа переходного температурного поля «Transient Thermal». 
При численном решении выполнялась калибровка параметров моделирования и граничных условий тепловой 
модели (таблица 1) для приближения модельных значений температуры к экспериментальным данным. 


Таблица 1 
Граничные условия (параметры) тепловой модели 


" Плотности тепловых Коэффициент 
Мощность тепловыделений| Тепловые потоки ов а 
Параметр a 


qvi, Вт/м? 42, Вт/м? Qi, Br Q2, Вт qni, Вт/м? qw, Вт/м? а, Br/(w?-?C) 


Значение 6 500 1 000 28 15 32 18 15 


В модуле «Ansys Mechanical» назначались начальные (начальная температура 790-24 °C) и граничные (таб- 
лица 1) условия для построенной сеточной модели, выполнялось построение температурного поля (рис. 2). При 


этом коэффициент теплопроводности принимался равным A=60,5 Вт/(м °С) и назначался таковым для конструк- 
ционной стали. 


Тепловая модель объекта включала два контактных соединения для электродвигателя и гильзы шпинделя с 


корпусом сверлильной головки и содержала 7 тепловых граничных условий. Построенная сеточная модель со- 
стояла из 16 309 элементов и 58 527 узлов. 


Общее время моделирования 21 600 секунд (6 часов) разбивалось на интервалы (1 час) и шаги At = 360 секунд 


(6 минут), всего № = 60 шагов, в пределах которых параметры тепловой модели (граничные условия) принима- 
лись как постоянные и независящие от времени. 


Моде! (В4) 
Hy Geometry 
#- „8 Materials 
Hy tk Coordinate Systems 
fyi} Connections 
yD Mesh 
=-- „М Transient Thermal (B5) 
vT- Initial Temperature 
via Analysis Settings 
B® Simplorer 
y®, Internal Heat Generation1 
y@, Internal Heat Generation2 
sa Heat Flow1 
sa Heat Flow2 
85 Heat Flux1 
85 Heat Fiux2 
Vf Convection 
S- Solution (B6) 
08-23 Solution Information 
s Temperature 
Bg «09 Simplorer 
/® ТетрРгоБе1 
8 TempProbe2 
i TempProbe3 
„® TempProbe4 


Рис. 2. Температурное поле объекта в «Ansys Mechanical» 


Расчет температурного поля требовался для генерации независимого пошагового отклика тепловой нагрузки, 
содержащего информацию об изменении температуры по времени (тепловые характеристики) для каждого пара- 
метра (таблица 1) тепловой модели в отдельности: 

УЕ (t) i=1,k, j=1,m, (5) 


где k — количество температурных датчиков; т — количество параметров тепловой модели. 
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Генерация пошагового отклика выполнялась в модуле «Ansys Mechanical» с применением расширения 
«Application Customization Toolkit (ACT)», которое поддерживает внедрение пользовательских сценариев. Это 
позволило разработать программный сценарий на языке программирования «Python» с целью автоматической 
генерации специального набора файлов, содержащих независимые пошаговые отклики [14]. 

Аналитическое решение выполнялось на основе метода пространства состояний путем построения модели 
тепловых характеристик по формуле: 

x= Ах+ Bu 


> 6 
y = Cx + Du (6) 


где x= ox / ôt — производная вектора х состояний по времени f; y и и — векторы выходных и входных данных 


соответственно; A, В, С, р — матрицы постоянных коэффициентов. 


В уравнении (6), вектор x-—(xi3»...xw) содержит переменные состояния, вектор входа 
u = (qui, 9, Оп, Q2, Чт, qm, Toi, Too, ..., Ток)! значения параметров тепловой модели (граничные и начальные усло- 
вия), вектор выхода y = (Ti, 75, ..., Ty)’ значения тепловых характеристик. 


Модель передаточной функции выражается следующим уравнением: 


т Е 
„(= У, His (uy abe, (7) 
где Н — матричная комплексная передаточная функция. 

Матричная передаточная функция Н получена путем применения к формуле (6) преобразования Лапласа, что 
выражается следующим уравнением: 

H(s)-C(sI- A4) B+D, (8) 
где 5 — комплексная переменная Лапласа; / — единичная диагональная матрица. 

Передаточная функция (8) отражает зависимость преобразования Лапласа выходной переменной 
Y(s) = H(s)U(s) от преобразования Лапласа входной переменной U(s) модели (6) при нулевых начальных 
условиях х(0) = xo = 0. При этом размерность матрицы передаточной функции H зависит от величины выхода 
k= 4 и размерности входной величины т = 10, что соответствует размерности исходного пошагового OT- 
клика (5). Поэтому модель (6) приводилась в соответствие с поведением исходной тепловой модели путем 
приближения ее передаточной функции (7) к пошаговому отклику (5) с применением метода векторной ап- 
проксимации [15]. 

Для построения матриц коэффициентов модели (6) выполнялся обратный переход от передаточной функ- 
ции (8) к модели в пространстве состояний. В этом случае, передаточная функция (8) принимает вид уравнения, 
в знаменателе которого содержится характеристический полином степени / = 4 (порядок системы), а в числителе 
полином степени Z = [- 1: 

z І 
G(s) = at P(s) = bg + du Q(s) = dg 2232077 (9) 
где a u b коэффициенты полиномов Q(s) и P(S) соответственно. 

Корни полиномов Q(s) и P(s) представляют собой полюса и нули передаточной функции (8) соответственно. 
Метод неопределённых множителей применялся к уравнению (9) для разложения каждого элемента матрицы Н 
на элементарные дроби. Обозначив полюса характеристического полинома через p;, получаем уравнение следу- 


G(s)=D+ Y". А 


i 
, 
5 — р; 


ющего вида: 


(10) 


2. Р(5) 
где R; E are 


Ранг матрицы R; обозначался через F; и выполнялось ее разложение на произведение двух матриц с полным 


— матрица размерностью (k X т); 4 — количество полюсов. 


рангом столбца и строки соответственно: 
kxr; pnxm ., = 
R; = C?" Bj". rank(R;) =". (11) 


Матрицы модели (6) являются диагональными, размерностью А", В", С", De" (п = 56) и содержат эле- 
менты, которые получены непосредственно из коэффициентов передаточной функции (10). 
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Матрицы A системы и В управления содержат следующие коэффициенты: 


ра 0 O° зв 0] bj 0 .. 0 
0 i 0 0 b e s. O 
P Pr ‚В| (12) 
0 0 рз 0 0 b, bj, 
[0 0 0 Pan | [0 O0 ud 
Матрицы С выхода и D прямой связи содержат следующие коэффициенты: 
Ca e Clg 0 sis 0 о... 0 0 .. 1.. 0 0 
C- alt m s at a 13 
0 es 0 (кл) C(k—1)g 0 ... 0 0 was Oa uu 1 0 ( ) 
0 .. 0 0 ses 0 Ck; e Chg 0. .. 0 .. 0 1 


Такой подход позволил выполнить переход от модели B форме передаточной функции (8) к модели B пространстве 
состояний (6). Алгоритмы перехода также реализованы в системе инженерных расчетов «MATLAB», коммерческой 
компании «The MathWorks» (США), в виде специальных функции «ss2tf()» и «tf2ss()». 

Для получения значений тепловых характеристик решалась задача Коши для системы обыкновенных диффе- 
ренциальных уравнений, поскольку в формуле (6) переменная X является производной от вектора состояний тем- 
пературного поля по времени t. Решение системы уравнений (6) было получено с помощью метода Рунге-Кутты 
четвертого порядка. 

Верификация построенной аналитической модели выполнялась путем проведения серии вычислительных экс- 
периментов с применением, разработанной в системе «MATLAB» прикладной программы, которая включала pe- 
ализацию метода Рунге-Кутты четвертого порядка (рис. 3). Вычислительные эксперименты проводились на пер- 
сональном компьютере (процессор AMD Ryzen 5 56000 with Radeon Graphics 2.30 GHz, оперативная память 
16,0 ГБ, тип системы 64-разрядная операционная система Windows 10 Pro версия 21H2), характеристики которого 
являются базовыми для современной вычислительной техники. 

Построенная аналитическая модель и применение метода Рунге-Кутты четвертого порядка для решения си- 
стемы дифференциальных уравнений позволило с высокой точностью вычислить значения тепловых характери- 
стик (рис. 3). Поскольку максимальные значения погрешности, то есть разность между модельными значениями 
тепловых характеристик, полученных с помощью численного (FE-Model) и аналитического (LTI-Model) реше- 
ний, не превысили 0,72 °С на всем интервале моделирования. 

Линейная и инвариантная по времени (Linear and Time-Invariant, LTI) компьютерная модель пониженного порядка 
(Reduced Order Model, ROM) разрабатывалась в модуле «Ansys Twin Builder» системы инженерного анализа «Ansys» 
на основе температурного поля объекта и сгенерированного в модуле «Ansys Mechanical» пошагового отклика. CCC 
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Рис. 3. Графики тепловых характеристик в системе «MATLAB» 
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Реализация цифровой модели тепловых характеристик на основе температурного поля с применением век- 
торной аппроксимации заключается в последовательном выполнении семи основных этапов: 

Этап 1. Импорт геометрической модели объекта и построение тепловой модели в системе инженерного анализа; 

Этап 2. Расчет температурного поля объекта на основе построенной тепловой модели; 

Этап 3. Генерация независимого пошагового отклика по результатам численного моделирования тепловых 
характеристик объекта; 

Этап 4. Применение алгоритма векторной аппроксимации для получения полюсов и нулей передаточной 
функции модели пространства состояний; 

Этап 5. Построение модели пространства состояний по известной передаточной функции; 

Этап 6. Разработка компьютерной модели тепловых характеристик; 

Этап 7. Реализация цифровой модели объекта. 

Таким образом предложенная цифровая модель, содержащая компьютерную модель для точной оценки теп- 
ловых характеристик сложных технических объектов проектирования, была получена путем последовательного 
выполнения всех представленных выше этапов. 

Результаты исследования. Разработанная компьютерная модель (Thermal SG400V_SML1) содержит 6 вхо- 
дов и 4 выхода (рис. 4), что позволяет выявить взаимосвязь между параметрами тепловой модели и значениями 
тепловых характеристик. В качестве входных данных компьютерной модели приняты объемная мощность теп- 
ловыделений (qui, 4), тепловые потоки (Qi, Q2) и плотности тепловых потоков (ni, 4,2). В качестве выходных 
данных выступают тепловые характеристики (T;—74). Компьютерная модель входит в состав цифровой MO- 
дели (рис. 4 а) тепловых характеристик, реализованной в модуле «Ansys Twin Builder». 

Значения параметров цифровой модели, представленные в табличном виде (рис. 4 в), подаются на вход ком- 
пьютерной модели (рис. 4 а) с помощью компонентов «STEP». Графический модуль демонстрирует значения 
тепловых характеристик (рис. 4 6), полученных на выходе компьютерной модели. 


Thermal SG400V SML 
RC 
34,0- 
T, 1 
34,074 
T; ] 
T, 2904 
Т, ] 
о аррар 
0 5 000 10 000 15 000 Lc 
a) 6) 
1 2 3 4 5 6 
Time, c 3 600 7 200 10 800 14 400 18 000 21 600 
qn2 18,0 18,0 18,0 18,0 18,0 18,0 
qi 6 500,0 6 500,0 6 500,0 6 500,0 6 500,0 6 500,0 
Q1 28,0 28,0 28,0 28,0 28,0 28,0 
Q2 15,0 15,0 15,0 15,0 15,0 15,0 
Qv2 1 000,0 1 000,0 1 000,0 1 000,0 1 000,0 1 000,0 
qnl 32,0 32,0 32,0 32,0 32,0 32,0 


Рис. 4. Реализация цифровой модели в системе «Ansys Twin Builder»: 
а — компьютерная модель тепловых характеристик; 6 — модуль графического представления результатов; 
в — модуль контроля входных данных компьютерной модели 


При разработке компьютерной модели задавались пределы размерности от 2 до 4 порядка, значения целевой 
ошибки £ = 5x10? и допуск для нулевого порядка £o = 2x10?. Остальные параметры устанавливались автомати- 
чески, поскольку метод векторной аппроксимации максимально автоматизирован в сравнении с остальными ме- 
тодами, которые поддерживаются в модуле «Ansys Twin Builder». 

В процессе разработки компьютерной модели модуль «Ansys Twin Builder» автоматически генерирует специ- 


альную матрицу ошибок аппроксимации М; = || у;()-у;(0 ||; во временной области, каждый элемент которой 
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отражает разность между значениями тепловых характеристик пошагового отклика (5) и передаточной функ- 
ции (7) модели в пространстве состояний, что выражается следующим уравнением: 
0,43 4,04 1,53 4,97 2,22 1,82 
_ 0,21 2,22 2,54 1,07 1,16 1,95 103 (14) 
2,66 1,63 0,39 0,26 3,43 1,82 
1,46 3,21 0,06 0,84 2,95 3,50 

В данном случае максимальная относительная ошибка не превышает значение £ = 4,97х103. При этом все 
остальные значения ошибок оказались меньше указанного 5 = 5x10? предела. Нулевое значение ошибки в мат- 
рице (14) означает, что вход был проигнорирован в связи с очень малым вкладом. 

Оценка точности расчета тепловых характеристик с применением предложенной цифровой модели выполня- 
лась путем сравнительного анализа результатов, полученных с помощью численного и аналитического решений. 
Сравнительный анализ проводился по критерию максимальной погрешности. Максимальная погрешность AT, 
то есть разность между значениями тепловых характеристик, полученных с помощью численного и аналитиче- 
ского решений для всех температурных датчиков, в каждый момент времени рассчитывается по формуле: 

AT nax = тахАТ,, (15) 


тах 
j-ln 


где АТ; = |Т,/- T;4 — погрешность; Т; ги 7; и — значения температуры конечно-элементной и цифровой моделей 


соответственно (j = 1, т); т — количество температурных датчиков. 

Для оценки точности цифровой модели выполнялся расчет погрешностей, результаты которого представлены 
в виде поверхности (рис. 5 а) и линейного графика (рис. 5 6) максимальной погрешности по времени. Поверх- 
ность (рис. 5 a) представляет собой расчетные значения погрешности AT; для каждого температурного датчика 
(ат, аТ2, dT3, dT4) в отдельности и в каждый момент времени на всем интервале моделирования. 
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Рис. 5. Погрешности расчета тепловых характеристик: а — погрешность для каждого температурного датчика; 
6 — максимальная погрешность для всех температурных датчиков 
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График (рис. 5 6) показывает, что для выбранного интервала времени максимальная погрешность AT, не 
превышает 0,1 °С. Погрешность расчета тепловых характеристик для каждого температурного датчика в отдель- 
ности (рис. 5 а) также не превышает заданных пределов. 

Обсуждение и заключение. Полученные результаты вычислительных экспериментов и проведенный срав- 
нительный анализ подтверждают эффективность применения предложенной цифровой модели, которая позво- 
ляет с высокой точностью AT. = 0,052 °С выполнить расчет тепловых характеристик, для дальнейшего прове- 
дения процедур анализа и идентификации тепловой модели. 

На температурное поле объекта проектирования влияет множество факторов, что осложняет определение номи- 
нальных значений тепловых граничных условий. Для решения данной проблемы в исследовании сначала анализи- 
руются ключевые технологии, задействованные в реализации цифрового двойника сложного технического объекта, 
затем строятся тепловая модель и компьютерная модель пониженного порядка LTI ROM переходного температур- 
ного поля. Разработанная компьютерная модель применяется в составе цифровой модели, которая позволяет полу- 
чить точную оценку тепловых характеристик объекта проектирования и тем самым повысить эффективность вы- 
полнения проектных процедур в процессе автоматизированного проектирования сложных технических объектов. 

Точность и эффективность вычислений компьютерной модели пониженного порядка и исходной тепловой 
модели полного порядка оценивается путем сравнительного анализа результатов моделирования. Результаты вы- 
числительных экспериментов показывают, что с точки зрения точности расчетов компьютерные модели пони- 
женного порядка и конечно-элементные модели полного порядка в основном сопоставимы по точности, а макси- 
мальная погрешность расчетов находится в пределах допустимого диапазона и не превышает 0,1 °С. 

Полученные в ходе вычислительных экспериментов результаты не противоречат результатам, представленным в 
источниках научной литературы по схожей тематике и позволяют сделать вывод о том, что применение предложенной 
цифровой модели эффективно для оценки тепловых характеристик сложных технических объектов в режиме реаль- 
ного времени, что является одним из важнейших условий для реализации технологии цифровых двойников. 

Однако изменение температурного поля сложного технического объекта по-прежнему зависит от многих фак- 
торов. Поэтому в дальнейших исследованиях необходимо разработать компьютерные модели тепловых дефор- 
маций и ввести эффективные алгоритмы оптимизации на основе искусственного интеллекта, чтобы обеспечить 
надежность результатов моделирования, полученных с помощью цифровых двойников. 
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Abstract 


Introduction. The dynamic loads during the start-up of a bridge crane can cause excessive stress in the structure and 
components, leading to potential safety hazards and increased wear and tear. To reduce the influence of the dynamic 
loads, various strategies can be implemented including optimization of the acceleration and deceleration profiles, using 
the soft start controls, implementing the vibration damping systems. It is vital to ensure that the proper crane maintenance 
and inspection protocols are in place. By reducing the impact of dynamic loads during the start-up, the overall 
performance and longevity of a bridge crane can be improved, ultimately enhancing safety and efficiency of the industrial 
operations. The present research offers a new approach to improving the efficiency and safety of industrial operations by 
providing a more precise account of the dynamic loads during the start-up of a bridge crane. The objective of this study 
is to develop a mathematical model for investigating the mechanical properties of the bridge cranes by analyzing the 
dynamic loads that occur during lifting operations. 

Materials and Methods. The development of the mathematical model was based on the kinetic model of the system, 
which included three connecting blocks and two flexible connections for a more accurate description of the bridge crane 
structure. Lagrange’s equations incorporating the information about the geometry and structure of a bridge crane were 
used. They made it possible to describe the motion of a system with the multiple elements and degrees of freedom. 
Processing and analysis of the results of the mathematical model were carried out in the MATLAB program using the 
Runge-Kutta method. 

Results. As a result of the research, a mathematical model was developed to study the dynamic loads affecting a bridge 
crane during lifting operations. Graphs describing the dependences of speed, acceleration, load, and rope angle over time, 
and their influence on the crane beam were plotted. The changes in these parameters over time, including their maximum 
values, were analyzed. The reasons for load changes and factors influencing the extension of lifting machines’ service 
life as well as reducing metal consumption during production thereof were identified. 

Discussion and Conclusion. The developed mathematical model and its numerical solution using the specialized software 
(MATLAB) allow for conducting the dynamic analysis of the bridge crane structures and determining the optimal design 
solutions. The analysis of the factors influencing the load changes leads to the conclusion that the use of this model can 
significantly reduce the load magnitudes and metal consumption, as well as increase the service life of lifting machines. 
The results obtained with the developed mathematical model and its numerical solution are useful for optimizing the crane 
structures, providing compliance with the operational requirements, and extending the service life of lifting machines. 
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Аннотация 

Введение. Динамические нагрузки во время запуска мостового крана могут вызывать избыточные напряжения в 
конструкции, приводя к потенциальным рискам и увеличению износа. Для снижения влияния динамических 
нагрузок можно применять различные стратегии, включая оптимизацию профилей ускорения и замедления, 
использование плавного пуска, внедрение систем амортизации. Важно обеспечивать исполнение правильных 
протоколов обслуживания и инспекции кранов. Путем снижения воздействия динамических нагрузок во время 
запуска можно улучшить общую производительность и долговечность мостового крана, повысив в конечном 
итоге безопасность и эффективность промышленных операций. Данное исследование предлагает новый подход 
к повышению эффективности и безопасности промышленных операций за счет более точного учета 
динамических нагрузок мостового крана при пуске. Цель работы — разработка математической модели для 
изучения механических свойств мостовых кранов путем анализа динамических нагрузок, возникающих во время 
подъемных операций. 

Материалы и методы. Разработка математической модели была выполнена на основе кинетической модели 
системы, включающей три соединительных блока и два гибких соединения для более точного описания 
конструкции мостового крана. Использованы уравнения Лагранжа, включающие информацию о геометрии и 
структуре мостового крана. Они позволили описать движение системы с несколькими элементами и несколькими 
степенями свободы. Обработка и анализ результатов математической модели были произведены в программе 
MATLAB с применением метода Рунге-Кутты. 

Результаты исследования. В результате исследования была разработана математическая модель для изучения 
динамических нагрузок на мостовой кран во время подъемных операций. Построены графики, описывающие 
зависимости скорости, ускорения, нагрузки и угла каната относительно времени и их влияние на балку крана. 
Проанализировано изменение этих параметров во времени, включая их максимальные значения. Определены 
причины изменений нагрузки и факторы, влияющие на увеличение срока службы и снижение металлоемкости 
при производстве подъемных машин. 

Обсуждение и заключение. Разработанная математическая модель и ее численное решение с использованием 
специализированного программного обеспечения (программа MATLAB) позволяют проводить динамический 
анализ конструкций мостового крана и определять оптимальные конструктивные решения. Анализ факторов, 
влияющих на изменение нагрузки, позволяет сделать вывод, что при использовании данной модели можно 
значительно снизить величину нагрузок и металлоемкость, а также увеличить срок службы подъемных машин. 
Результаты, полученные при помощи разработанной математической модели, и ее численное решение полезны 
при оптимизации конструкции кранов, обеспечении соответствия операционных требований и продлении срока 
службы подъемных машин. 


Ключевые слова: мостовой кран, динамическая нагрузка, кинетическая модель, подъем грузов, динамический анализ 
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1653-2024-24-2-190-197 


Introduction. The analysis of dynamic processes in the mechanical part plays an important role in the development 
of new overhead cranes and modernization of existing ones in order to reduce loads on control devices and extend their 
service life [1]. 

The bridge crane is subjected to dynamic loads during non-static operations, such as acceleration and braking. Analysis 
of these processes allows identifying hidden impacts on the dynamic behavior of the bridge crane. Therefore, it is 
paramount for the researcher to make an optimal design choice to reduce these loads, ensuring that the crane can meet the 
required operating conditions [2]. 

In [3], a dynamic model of a crane lifting system was developed, using which an accurate direct numerical integration 
method was proposed for calculating the dynamic loads of the system. 
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Papers [4, 5] studied dynamic loads in a metal structure, taking into account fatigue of the metal material. However, 
the researcher neglected the impact of forces from the drive of the lifting mechanism operating with artificial parameters. 

In [6], dynamic loads in a bridge crane were determined during the operation of a lifting mechanism when hoisting a 
load suspended on a rope. The most important case studied was the effect of dynamic loads on the crane when removing 
a load from a solid foundation, at the moment of lifting-off. 

Steel structures of bridge cranes experience non-stationary loads with different stress amplitudes and asymmetry of 
the working cycle [7]. To study the real load of bridge cranes under typical operating conditions, constant recording of 
their stress state 1s required, which 1$ labor-intensive [8]. Therefore, statistical processing of the obtained results 1$ used 
to assess the loading elements of metal structures of bridge cranes [9]. This involves changing individual components of 
the total load of metal elements, such as the gravity of the load being lifted, the angle of rotation of the load, and weather 
loads, and then summing them according to the laws of probability theory [9,10]. This approach is less labor-intensive 
than a comprehensive study of loading under typical operating conditions. However, determining the probabilistic 
characteristics of individual random loads also takes time, so the method for calculating load combinations 1s widely used 
in crane construction [11]. 

In [12], it was found that during the stage of selecting the rope slack, the value of the stator current of the electric 
motor of the lifting mechanism did not depend on the mass of the suspended load. However, as the load increased, the 
time of its rise also increased, and at the stage of separating the load from the surface, the amplitude values of the current 
increased. Furthermore, a noticeable difference appeared after five periods of mains voltage from the beginning of the 
stage. However, the researcher neglected the significant influence of forces arising from the drive of the lifting mechanism 
working with artificial elements. 

By using the developed mathematical model, the research aims at optimizing the design of bridge cranes through 
studying the dynamic loads that occur during lifting operations, ensuring that the crane has the ability to meet the required 
operating conditions. 

Materials and Methods 

Kinetic Model of the System under Study. When developing a kinetic model of the system under study, it can be 
represented that the construction of a bridge crane for a given motion form consists of three connecting blocks and two 
flexible joints, and has the form shown in Figure 1: 


Fig. 1. Bridge crane kinetic model 


тв — mass of the main bridge, moving along the X-axis; тр — mass of the lifting mechanism, moving along the Y-axis; 
тс — payload; Кв and Св — elasticity modulus and damping coefficient of the main bridge; Kr and Cr — elasticity 
modulus and damping coefficient of the ropes; a — swing angle of lifting mechanism winding. 

Derivation of Kinetic Equations Reflecting the Motion of the Dynamic Model 

Mathematical equations representing the motion of the dynamic model are derived from the partial differential 
Lagrange equation, which is considered to be one of the best methods used specifically in cases where the system consists 
of more than one element and when it has several degrees of freedom. 
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where T — kinetic energy; U — potential energy; D — damping energy; О — external forces affecting the whole system; 
qi — common system of coordinates; i — degrees of freedom of the model under study. 
Kinetic energy equation: 


T =0.5тв22 +0.5mpz% --0.5mg 22, --0.5J -à?. (4) 
Potential energy equation: 

U =0.5K ,z2 +0.5К, (za -ra za). (5) 

Damping equation: 

2 

D = 0.5C525 +0.5С, (2$ -r& € 25) . (6) 

System equations: 
(ть * ть)» * Kgzg - K,(zg ^ra * zo) Cg25 * Cr(Zg -rà& +26) =0, (7) 
"mois + К, (во такав) +С, оға, - (8) 
'J-&- K,r(a—zg -zg) +С, (0-26 -259) Ma -М,, (9) 


where Ma, M. — resistance momentum and lifting mechanism momentum; J — inertia of rotating mass of the lifting mechanism. 

Numerical Solution of a Mathematical Model Using the Fourth-Degree Equation of Runge-Kutta 

The numerical solution of equations (7-9) was obtained using the Runge-Kutta method in the MATLAB program. 
The dynamic equations were derived within the program, incorporating the input data and a set of commands to process 
these equations. The resulting graphs illustrate the interconnections between the various blocks and components of the 
crane structure under consideration. 

The Conditioning of the Input Values Required to Solve the Model 

The study was conducted on an ACE type bridge crane consisting of three parts (Fig. 2): 

— lifting trolley; 

— main bridge, which supports the lifting mechanism; 

— end trucks, which support the main bridge. 


Fig. 2. Main parts of a bridge crane ACE 


To align the operation of the bridge crane with a standard set of coordinate axes, the following assumption was made: 

— the lifting trolley functions as a unit responsible for raising a load along the Z-axis and allows for horizontal motion 
along the Y-axis relative to the main bridge, which 1$ the axis along which the crane moves. 

A study was conducted on a prototype bridge crane with a lifting capacity of 10 tons and a width of 21.5 meters. The 
following characteristics were considered: 

— payload: тс = 10,000 kg; 

— mass of the two main bridges: mg = 8,100 kg; 

— mass of the lifting mechanism with the trolley: mp — 700 kg. 

When determining the input values, all the laws of designing the structures of lifting devices were followed, and the 
connections of all components were taken into account. The main factors considered were: 

— power of the drive of the lifting mechanism along the Z-axis; 

— stiffness coefficient of the steel structure (Кв); 

— stiffness coefficient of the rope (К,); 

— damping coefficient of the metal frame (Св); 

— damping coefficient of the ropes (С,). 
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Table 1 
Shows the elements of the design structure of the bridge crane that was studied. 

Element name Value Unit of measurement Notation 
Load 10 Ton тс 
Crane mass 8,100 Kg тв 
Trolley mass 700 Kg mp 
Crane length 21.5 m L 
Height of lift 5 m H 
Beam device degree 2 e A 
Gear box ratio 4.5 — im 
Coil radius 0.25 m R 
Rope diameter 16.5 mm di 
Engine power 30 kW Nn 
Speed of the engine rotor core 905 r.p.m Nn 
Rope stiffness coefficient 11 169.8 N/mm Кв 
Rope damping coefficient 23 934.4 N/mm К, 
Rope damping coefficient 83.37 N.sec/m С, 
Metal frame damping coefficient 30.29 N.sec/m Св 


Research Results 


Study of the Model Operation Using a Computer 
Structural Behavior of a Metal Bridge Crane under Momenta Loads 
The structural behavior of a metal bridge crane during the process of lifting a load is depicted by three curves (Fig. 3). 
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Fig. 3. Coordinates of vertical motion, speed and acceleration of the COG of the bridge cranes versus time curves 


In Figure 3, the first curve shows the coordinates of the bridge crane along the ordinate axis as a function of time. 
During lifting, some vibration of the crane's metal structure is observed for 5—10 seconds, which then stabilizes and does 


not affect its rigidity. 


The second curve in this figure depicts the change in the speed of the center of gravity of the bridge crane over time. 
During lifting the load, the speed initially increases and then gradually decreases until stabilization. This indicates that 
the center of gravity of the bridge crane assumes a stable position within the specified time period, during which the 


vibration stabilizes. 


The third curve reflects the change in the center of gravity of the bridge crane over time. It is noteworthy that the 
acceleration value at the moment of lifting the load is 0.07 m/s?, with the dynamic load reaching its maximum.? 


Antypas IR. Modeling the Dynamic Loads Affecting a Bridge Crane during Start-Up 


Load Curves 
Figure 4 shows three curves that reflect motion of the load during the crane operation at the time of lifting. 
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Fig. 4. Coordinates of the speed and acceleration of the load lifting versus time curves 


In Figure 4, the first curve reflects the change in the position of the center of gravity of the load over time during its 
lifting to a certain height, calculated by the program, by rotating the drum by 180 degrees and then stopping. The height 
of the load suspension indicated on the graph 1$ 0.78 m. 

The second curve on this graph shows the change in the lifting speed of the load over time. When the drum rotates, 
the lifting speed of the load initially increases, then gradually decreases until it stabilizes. 

The third curve represents the change ш the acceleration of the load lifting over time. At the moment of lifting the 
load, the acceleration reaches a maximum value of 0.26 m/s? in 0.8 seconds, and then stabilizes. 

From the analysis of these three curves, it can be noticed that the stabilization time of the crane operation is almost constant. 

Coil Angle Curve 

Figure 5 shows the dependence of the winding rotation on time in degrees. 
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Fig. 5. Coil Angle Curve 


Figure 5 shows the change in the angle of rotation of the drum over time. The rotation angle stabilizes when reaching 
a value of 180 degrees, after which it remains stable, meaning it repeats. 
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Discussion and Conclusion. Analysis of the above graphs leads to the following conclusions. 

The mathematical model and algorithms allow for a detailed study of the motion of a bridge crane at all stages. 
Adjusting the winch operation and improving the metal structure of the crane have many positive aspects, including 
reducing dynamic displacements in the metal frame of the crane and transmission elements (such as clutch, gearbox, 
motor, and pulleys), as well as in the hoisting ropes. This reduction in dynamic loads leads to a decrease in rapid wear of 
these elements. Additionally, cost savings on maintenance and the development of an optimal metal structure design are 
achieved, as well as an increase in the crane's service life. This is evident from changes in the center of gravity of the 
load, speed, and acceleration during its lifting over time. It has been established that the height of the load during lifting, 
which is directly related to the length of the hoisting ropes, as well as the mass of the lifted load and the design of the 
main bridge, including its shape and dimensions, significantly influence the dynamic loads experienced by the structure. 

The mathematical model offers a new approach to improving the efficiency and safety of industrial operations by 
providing a more precise understanding and accounting for dynamic loads during the start-up of a bridge crane. 
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Аннотация 

Введение. Безопасность судоходства и разработок подводных месторождений полезных ископаемых требуют 
точного обнаружения различных подводных объектов. В литературе рассматриваются вопросы отслеживания их 
перемещений и траектории движения. Предлагаются методы гидролокации, обеспечивающие высокую точность 
позиционирования подводных объектов. Отмечена высокая точность пеленга стереодатчиков с ультракороткой 
базой. Однако такое оборудование чувствительно к частоте дискретизации сигналов, что вызывает «шум дискре- 
тизации». В открытом доступе нет публикаций, посвященных решению этой проблемы. Представленное иссле- 
дование призвано восполнить данный пробел. Цель работы — изучение возможности получения данных, уточ- 
няющих информацию о пеленге подводных объектов за счет использования фазовой информации отраженных 
зондирующих сигналов и дополнительной процедуры передискретизации исходных данных. 

Материалы и методы. Местоположение объекта определяли с помощью экспериментального комплекса для 
исследования гидроакустических датчиков, созданного В.А. Широковым и В.Н. Милич в Удмуртском федераль- 
ном исследовательском центре Уральского отделения Российской академии наук. Использовали стереодатчик с 
малой базой (30 мм) по сравнению с расстоянием до объекта (=800-900 мм). Для обработки данных применяли 
методы цифровой фильтрации и математический аппарат корреляционного анализа отраженных гидроакустиче- 
ских сигналов, полученных фазовым методом. 

Результаты исследования. Представлены итоги сопоставления двух способов определения пеленга на объект: 
по разности времени прихода передних фронтов импульсов и по максимуму кросс-корреляционной функции 
(ККФ). Графически показано изменение пеленга при движении объекта. Использование переднего фронта сиг- 
нала обусловило небольшие выбросы значений вдоль всей кривой пеленга (менее 0,12 рад). При максимуме ККФ 
выбросы фиксировались лишь в некоторых областях, но были довольно значительными (около 0,17 рад). Пока- 
зано, как выбрать точки, соответствующие более гладкой и валидной траектории объекта, и как работать с оши- 
бочными точками. Представленный метод устранения ошибки можно реализовать программно. При квазигармо- 
ничном сигнале редкие измерения исходного сигнала интерполируются частыми вычисленными значениями. 
Благодаря такому виртуальному увеличению частоты дискретизации (передискретизации) можно фиксировать 
промежуточные показатели в оцифрованных исходных данных. Интерполяция значений сигнала кубическим 
сплайном позволила получить 20 точек на 1 период сигнала вместо 5 точек в исходном варианте. В этом случае 
более корректна траектория, сформированная с максимумом ККФ. 

Обсуждение и заключение. Задачу пеленгации можно решить с точностью, необходимой для практического 
применения. Учет фактора гладкости и непрерывности траектории движения объекта позволяет качественно кор- 
ректировать выбор максимума кросс-корреляционной функции сигналов стереодатчика. Предложенные методы 
обладают большим потенциалом для разработки систем подводного видения. 
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Ключевые слова: определение местоположения объекта в гидросреде, пеленгация фазовым методом, передние 
фронты импульсов, кросс-корреляционная функция, шум дискретизации 
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Abstract 

Introduction. Safety of navigation and development of underwater mineral deposits require the accurate detection of 
various underwater objects. The literature discusses the issues of tracking their motion and trajectory. Sonar methods are 
proposed to maintain high accuracy of underwater object positioning. High accuracy of the bearing of stereo sensors with 
an ultrashort base is noted. However, this equipment 1s sensitive to the sampling rate of the signals, which causes 
"sampling noise". There are no publicly available publications dedicated to the solution to this problem. The presented 
study is designed to fill this gap. This work is aimed to study the possibility of obtaining data clarifying information about 
the bearing of underwater objects through using the phase information about echoed probing signals and an additional 
procedure for resampling the source data. 

Materials and Methods. The location of the object was determined using the experimental complex for studying 
hydroacoustic sensors created by V.A. Shirokov and V.N. Milich at the Udmurt Federal Research Center, the Ural Branch 
of the Russian Academy of Sciences. A stereo sensor with a small base (30 mm) was used compared to the distance to 
the object («800—900 mm). Digital filtering methods and mathematical apparatus of correlation analysis of return 
hydroacoustic signals obtained by the phase method were used for data processing. 

Results. The results of comparing two methods for determining the bearing on an object are presented: by the difference 
in the time of arrival of the pulse-leading edges and by the maximum of the cross-correlation function (CCF). The change 
in bearing as the object moves, is graphically shown. The use of the leading edge of the signal caused small outliers of 
values along the entire bearing curve (less than 0.12 rad). At the maximum CCF, emissions were recorded only in some 
areas, but they were quite significant (about 0.17 rad). It showed how to select points corresponding to a smoother and 
more valid object trajectory, and how to work with erroneous points. The presented method of error correction can be 
implemented programmatically. With a quasi-harmonious signal, rare measurements of the original signal are interpolated 
by frequent calculated values. Thanks to this virtual increase in the sampling rate (oversampling), intermediate indicators 
can be recorded in the digitized source data. Interpolation of the signal values by a cubic spline allowed us to obtain 20 
points for 1 period of the signal instead of 5 points in the original version. In this case, the trajectory formed with the 
maximum CCF is more correct. 

Discussion and Conclusion. The direction-finding problem can be solved with the accuracy required for practical 
application. Taking into account the factor of smoothness and continuity of the object's trajectory makes it possible to 
qualitatively correct the selection of the maximum of the cross-correlation function of the stereo sensor signals. The 
proposed methods have great potential for the development of underwater vision systems. 


Keywords: determination of the location of an object in a hydraulic environment, phase direction finding, pulse-leading 
edges, cross-correlation function, sampling noise 
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Введение. Обеспечение безопасности судоходства и работы на подводных месторождениях полезных иско- 
паемых требуют качественного обнаружения подводных объектов и отслеживания их перемещений [1]. В систе- 
мах подводного наблюдения применяются гидроакустические датчики. Они улавливают сигнал, отраженный от 
объекта, и позволяют вычислить его местоположение методом трилатерации [2]. В этом процессе каждый датчик 
предоставляет информацию об интервале времени прохождения зондирующего сигнала, отраженного от объекта. 
В случае использования нескольких датчиков [3] возникает возможность решения задачи обратной простран- 
ственной засечки и определения координат наблюдаемого объекта [4]. Для повышения точности измерений пред- 
ставляется перспективным использование в качестве приемника стереодатчиков с ультракороткой базой [5], поз- 
воляющих получать фазовую информацию [6] и определять пеленг на объект [7]. Изучены возможности исполь- 
зовании гидролокационных данных об удаленных подводных целях, а также о пеленге подводных объектов [8] 
за счет использования фазовой [9] или частотной [10] информации отраженных зондирующих сигналов. 

При этом известно, что процедура дискретизации сигнала обусловливает ошибки при определении парамет- 
ров траектории движения подводных объектов [11]. Это относится как к дальномерным измерениям, так и к про- 
цедурам определения пеленга. Оцифровка аналогового сигнала датчика, необходимая для дальнейшей цифровой 
обработки, вносит так называемую «погрешность дискретизации». Погрешность, вызванная квантованием ам- 
плитуды сигнала, оценивается количеством уровней квантования аналого-цифрового преобразователя. Погреш- 
ность, вызванная дискретизацией по времени, пропорциональна величине интервала времени квантования и ско- 
рости изменения напряжения. Поэтому частоту дискретизации сигнала стремятся сделать как можно выше верх- 
ней частоты самого сигнала. Однако увеличение частоты дискретизации во многих случаях ограничивается воз- 
можностями регистрирующей аппаратуры: процессора, аналого-цифрового преобразователя, каналов передачи 
данных, устройствами хранения данных. Ограниченность возможностей дискретизации во времени при выпол- 
нении траекторных измерений не позволяет точно фиксировать экстремальные значения траектории (дальность 
и пеленг). Как следствие, снижается точность результатов измерений, появляются выбросы на фиксируемых тра- 
екториях. Для устранения этого недостатка необходимо изучить потенциал передискретизации оцифрованного 
гидроакустического сигнала. В открытом доступе нет публикаций, посвященных решению данной проблемы. 
Имеющийся пробел восполняют материалы этой статьи. 

Цель представленного исследования состоит в оценке возможностей определения пеленга на объект с исполь- 
зованием фазовой информации и в отработке фазового метода пеленгования с применением передискретизации 
оцифрованного гидроакустического сигнала. 

Материалы и методы. Рассматривается алгоритм определения координат с применением пеленгации фазо- 
вым методом. Приводятся результаты его тестирования для объекта, движущегося по круговой траектории в экс- 
периментальном бассейне. 

Постановка эксперимента. При проведении опытов для определения местоположения объекта в гидросреде 
при пеленгации фазовым методом задействовали экспериментальную установку, созданную В.А. Широковым и 
В.Н. Милич в Удмуртском федеральном исследовательском центре Уральского отделения Российской академии 
наук [12]. Это лабораторный измерительный комплекс с линейной аквасредой в виде протяженного цилиндриче- 
ского резервуара (гидроволновода) и бассейном, оснащенным системой генерации испытательных гидроакусти- 
ческих сигналов. Ниже детально описаны элементы установки, которые использовались в опытах. 

1. Излучатель гидроакустического сигнала с известными координатами 5 (рис. 1). В работе используется ам- 
плитудно-модулированный сигнал [13]. 


M 


Рис. 1. Расположение принимающих датчиков (стереодатчик с двумя приемниками 51 и 512) и излучателя (5) B плоскости 
экспериментального бассейна (xy). Здесь 4 — база стереодатчика; D(O; st) — расстояние между стереодатчиком и объектом; 
а — пеленг на объект; плоскость (x у’) образована плоскостью стереодатчика и нормалью и к ней 
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2. Отражающий объект. В эксперименте в качестве тестового используется протяженный цилиндрический 
объект — медный провод, диаметр которого (0,15 мм) значительно меньше длины акустической волны (1,5 мм). 
В представленной работе задача решается на плоскости, поэтому для составления вычислительных процедур 
расчета координат данный объект можно считать точечным. 

3. Стереодатчик, преобразующий гидроакустический сигнал в электрический. Известны координаты каждого 
из двух приемников стереодатчика (511, 512 на рис. 1). 

4. Аппаратно-программный комплекс, выполняющий усиление, оцифровку и обработку сигналов датчиков. 

Результатом измерений является файл регистрации принимаемого двухканального сигнала. 

Алгоритм обработки результатов измерений. Использование стереодатчика с малой базой (30 мм) по срав- 
нению с расстоянием до объекта (= 800—900 мм) позволяет применить упрощенные формулы. 

Итак, рассчитаем координаты исследуемого объекта в системе x'y', образованной плоскостью стереодатчика 
и нормалью к ней. Для этого задействуем алгоритм, включающий пять этапов. 

1. Предварительная обработка: 

— удаление слепой зоны (первые М отсчетов); 

— фильтрация сигнала. 

В работе применяется высокочастотная фильтрация. 

2. Выделение из всего сигнала информативных фрагментов impl, imp2 (в двух каналах), содержащих им- 
пульсы, отраженные от исследуемого объекта. Для этого выполняется пороговая обработка по амплитуде А сиг- 
нала. Порог принимается равным рхтах(А) (р = 0,5). 

Такая обработка позволяет определить передний фронт импульса (BeginIndex), отраженного от объекта. Из двух 
каналов сигнала вырезаются участки [BeginIndex-0,5xLenSignal, BeginIndex+LenSignal]. Здесь LenSignal — длина car- 
нала, подаваемого на вход. В результате получаем два сигнала impl и ітр2 длиной 1,5xLenSignal. В каждом из них 
содержится импульс, отраженный от объекта. 

3. Расчет расстояния D(O; st) между объектом и стереодатчиком. 

— Определение расстояния, пройденного импульсом от излучателя до объекта и от объекта до датчика при- 
ema, по формуле D(S,O,st)=BeginIndexxdtxC. Здесь C=1475 м/с — скорость распространения акустической волны, 
dt = 0,2x10—6 c — интервал дисктретизации. 

— Расчет значения: 

D(S,O,st)—D(S, st) 


D(O, st) = 3 (1) 
4. Определение пеленга а на объект [14]: 
a = arcsin Ат (2) 


где À — длина волны; F — частота сигнала; dt — интервал дискретизации; т — разница прихода отраженного 
импульса Ha два стереодатчика (в отсчетах); d — база стереодатчика (в представленном исследовании A = 1,5 мм, 
Е= 1 МГЦ, dt = 0,2 мкс, d = 30 мм). 

5. Вычисление координат анализируемого объекта в системе, образованной плоскостью стереодатчика и его 
нормалью (рис. 1), по формулам: x —D(O, st) *cos(a); y =D(O, sf) xsin(a). 

Определение пеленга на объект. Примем условие малости стереобазы по сравнению с расстоянием от дат- 
чика до объекта. В этом случае стереодатчик позволяет определить пеленг на объект с помощью информации о 
разности фаз принятых импульсов, ранее отраженных от объекта, а также по разности времени прихода передних 
фронтов импульсов, поступивших на стереодатчик (рис. 2). 


AB, 


Th\ 
Th2 


Рис. 2. Определение разности моментов прихода ds стереосигналов, которые обозначены сплошной imp 1 (f) 
и пунктирной imp2(t) линиями. Thl, Th2 — пороги для первого и второго сигнала стереопары соответственно; 
dt — интервал дискретизации; А — значение сигнала B вольтах; f — время 
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Для реализации такого подхода выполняются описанные ниже процедуры. 

1. Идентификация моментов прихода передних фронтов стереосигналов с помощью пороговой обработки [15] 
с порогом рхтах(А). В данной работе коэффициент р принят равным 0,5 [16]. 

2. Вычисление разности ds (количество интервалов df) между передними фронтами импульсов, зарегистриро- 
ванных в двух приемниках стереодатчика. 

3. Вычисление пеленга a по формуле (1). За т принимаем значение ds. 

При таком способе вычисления сдвига фазы передний фронт сигнала будет главным индикатором импульса, 
отраженного от объекта. Не учитывается внутренняя структура сигнала, которая может иметь решающее значе- 
ние для корректного определения пеленга на объект. Учтем схожесть принимаемых сигналов и рассмотрим в 
работе другой подход к вычислению пеленга — на основе кросс-корреляционной функции (ККФ) стереосигна- 
лов. Логично определять сдвиг фазы по максимуму ККФ. Ниже описаны элементы алгоритма вычисления пе- 
ленга с помощью ККФ. 

1. Расчет ККФ импульсов, поступивших на два приемника одного стереодатчика: г = xcorr(impl, imp2). Зна- 
чение функции xcorr при смещении m второго сигнала относительно первого вычисляется следующим образом: 
N=m-1 

impl,,,,x imp2,, т>0, 
XCOIT;,51 imp2 (m)= =, 
ХСОГТ2 тр (7), m < 0. 

Здесь № — длина impl и imp2, -N«m«N. 

2. Определение максимума ККФ и соответствующего ему сдвига dc (в orcuerax). 

3. Вычисление пеленга а по формуле (1) с учетом, что за т принимается величина dc. 

Результаты исследования 

Тестирование методов определения пеленга. Представленные способы определения пеленга на объект (по 
разности времени прихода передних фронтов импульсов и по максимуму ККФ) протестировали для объекта, 
движущегося по круговой траектории в экспериментальном бассейне. Использовался гармонический сигнал дли- 
тельностью 7 периодов и частотой 1 МГц. Оцифрованные стереосигналы регистрировались в 256 точках траек- 
тории движения объекта. Частота оцифровки — 5 МГц. 

На рис. 3 показано изменение значений пеленга по мере движения объекта по окружности для двух предло- 
женных подходов (разница передних фронтов, максимум ККФ). 
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Рис. 3. Изменение пеленга а на объект, совершающий круговое движение. Показатели получены при использовании двух 
подходов (красный маркер — разница времени прихода передних фронтов стереосигналов, синий маркер — максимум 
ККФ). MN — номер измерения, а — пеленг на объект 
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При использовании переднего фронта сигнала (рис. 3, красный маркер) почти вдоль всей кривой изменения 
пеленга на объект наблюдаются небольшие выбросы значений. При использовании максимума ККФ (рис. 3, си- 
ний маркер) только в некоторых областях присутствуют выбросы значений пеленга, HO они достаточно большие. 
Возможно, это объясняется особенностями движения объекта в области, где зондирующий и отраженный сиг- 
налы пересекают зону турбулентности, вызванной движением объекта. При нахождении фазового сдвига стерео- 
сигналов с помощью максимума ККФ возникают дополнительные ошибки, что затем отражается в вычисленной 
траектории объекта. Также выбросы в вычисленных значениях пеленга могут быть обусловлены недостаточно 
высокой частотой дискретизации. Значения этих ошибок существенно выше показателей выбросов, наблюдае- 
мых для метода с использованием переднего фронта импульса. 

На рис. 4 а представлены траектории движения объекта, вычисленные с использованием глобального макси- 
мума ККФ и его соседних локальных максимумов. 
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Рис. 4. Анализ траекторий кругового движения объекта: а — траектории, вычисленные с использованием глобального 
максимума ККФ (черный маркер), локального ближайшего правого максимума ККФ (оранжевый маркер), локального 
ближайшего левого максимума ККФ (зеленый маркер); 6 — пример ККФ, глобальный максимум которой соответствует 
корректному значению координат объекта; в — пример ККФ, в которой корректному значению координат объекта 
соответствует левый локальный максимум ККФ 


Используя априорную информацию о плавности перемещения экспериментального объекта, можно выбрать 
точки, которые соответствуют более гладкой и валидной траектории. Такой подход применяется в других рабо- 
тах [17]. На рис. 4 а при x’e[800;870], y >0 корректные значения находятся на кривой для траектории, рассчитан- 
ной по левому локальному максимуму ККФ. Есть также ошибочные точки. Если их заменить точками траекто- 
рии, полученной в результате использования правого локального максимума ККФ, можно вычислить корректное 
значение точки траектории (x 800, у’<0 на рис. 4 а). Данный графический метод устранения ошибки может быть 
реализован программно. 

Таким образом, в некоторых точках траектории корректный сдвиг фазы стереосигналов, необходимый для 
вычисления пеленга, может соответствовать не глобальному максимуму ККФ, а одному из соседних локальных 
экстремумов (рис. 4 6, в). Такое смещение максимумов ККФ может быть обусловлено недостаточно высокой ча- 
стотой дискретизации сигналов. Этим же объясняется ступенчатый характер результирующей траектории. 

С учетом квазигармоничности сигнала предлагается использовать интерполяцию редких измерений исход- 
ного сигнала частыми вычисленными значениями. Такое виртуальное увеличение частоты дискретизации 
(передискретизации) дает возможность зафиксировать промежуточные значения в оцифрованных исходных дан- 
ных. Интерполяция значений сигнала кубическим сплайном позволила получить 20 точек на 1 период сигнала 
вместо 5 точек в исходном варианте. На рис. 5 приведены расчетные траектории движения объекта, сформиро- 
ванные с предварительной передискретизацией сигналов. 
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Рис. 5. Расчетные траектории движения объекта, полученные для стереосигналов с виртуально повышенной частотой 
дискретизации: синий маркер — используется разница времени прихода передних фронтов сигналов; 
красный маркер — используется максимум ККФ) 
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Красная линия — это траектория движения объекта, полученная с использованием максимума ККФ для pac- 
чета пеленга по предварительно интерполированным данным. Она представляется более корректной по сравне- 
нию с аналогичной траекторией, полученной с использованием разницы передних фронтов отраженных импуль- 
сов (синяя линия). 

Обсуждение и заключение. Анализ особенностей определения пеленга на подводный объект с использова- 
нием взаимной фазовой информации сигналов дифференциального стереодатчика дает возможность сделать сле- 
дующие выводы: 

— разность времени прихода передних фронтов импульсов, зарегистрированных стереодатчиком, позволяет получать 
значения пеленга, показывающие направление на объект. Качество результата зависит от принятой величины порога (ис- 
пользуемого при определении переднего фронта импульса) и изменчивости амплитуд принимаемых сигналов; 

— применение подхода с использованием максимума кросс-корреляционной функции при недостаточной ча- 
стоте дискретизации исходных сигналов приводит к существенным выбросам в расчетной траектории (рис. 3); 

— передискретизация сигналов с помощью интерполяции исходных квазигармонических сигналов с исполь- 
зованием максимума ККФ для расчета пеленга и координат объекта позволяет получить гладкую траекторию 
(рис. 5), соответствующую плавному перемещению объекта по кругу. Такой подход требует значительно боль- 
ших вычислительных затрат в сравнении с другими, рассмотренными в этой статье. Это может повлиять на ско- 
рость обработки сигналов в режиме реального времени. 

Таким образом, исследованные возможности уточнения фазового метода пеленгации (передискретизация 
оцифрованных сигналов и использование максимума ККФ) позволяют эффективно оценивать траекторию дви- 
жения подводного объекта. 

Описаны условия получения гладкой корректной траектории движения отслеживаемого объекта. Этой цели 
служит передискретизация оцифрованных с недостаточной частотой дискретизации интерполированных квази- 
гармонических стереосигналов. Вместе с этим задействуется метод расчета пеленга, использующий максимум 
кросс-корреляционной функции сигналов стереодатчика. 

Данные на рис. 5 подтверждают, что задачу пеленгации можно решить с точностью, необходимой для прак- 
тического применения. Учет фактора гладкости и непрерывности траектории движения объекта позволяет одно- 
значно корректировать выбор максимума кросс-корреляционной функции сигналов стереодатчика (рис. 4). 

Предложенные методы имеют большое значение для разработки систем подводного видения. 
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